- Email: [email protected]
Functional modeling of neural–glial interaction
BioSystems 89 (2007) 84–91
Functional modeling of neural–glial interaction D.E. Postnova , L.S. Ryazanovaa , O.V. Sosnovtsevab,∗ a
Physics Department, Saratov State University, Astrakhanskaya Street 83, Saratov 410026, Russia Department of Physics, The Technical University of Denmark, 2800 Kongens Lyngby, Denmark
b
Received 24 December 2005; accepted 15 April 2006
Abstract We propose a generalized mathematical model for a small neural–glial ensemble. The model incorporates subunits of the tripartite synapse that includes a presynaptic neuron, the synaptic terminal itself, a postsynaptic neuron, and a glial cell. The glial cell is assumed to be activated via two different pathways: (i) the fast increase of intercellular [K+ ] produced by the spiking activity of the postsynaptic neuron, and (ii) the slow production of a mediator triggered by the synaptic activity. Our model predicts the long-term potentiation of the postsynaptic neuron as well as various [Ca2+ ] transients in response to the activation of different pathways. © 2006 Elsevier Ireland Ltd. All rights reserved. Keywords: Glion; Neural–glial interaction; Functional modeling
1. Introduction During the last decade, application of a variety of new tools has allowed us to uncover the diverse dynamic roles of glial cells. Glial cells are active partners to neurons in processing information and synaptic integration (Field and Steven-Graham, 2002; Bonvento et al., 2002). Thus, studying the interaction between the neurons and glial cells has become an important problem in neurophysiology and cell biophysics. The modeling of neural–glial interaction may help to understand the transition between normal and different pathological states. An adequate mathematical model of neural–glial interaction should account for all relevant ionic currents, processes of production, propagation of neuromediators, etc. Such a model inevitably will be high-dimensional with a large number of control parameters. While the
∗
Corresponding author. Tel.: +45 4525 3106; fax: +45 4593 1669. E-mail addresses: [email protected] (D.E. Postnov), [email protected] (L.S. Ryazanova), [email protected] (O.V. Sosnovtseva). URL: http://www.fys.dtu.dk/∼olga.
main pathways of neural–glial interactions seem to be quite common, the details can vary from species to species. Hence, there will be a number of speciesoriented mathematical models. Nadkarni and Jung (2003, 2004) first introduced the concept of a “dressed neuron”, representing a set of equations that describe the activation of an astrocyte by a neuron via the inositol 1,4,5-triphosphate (IP3 ) production that triggers the calcium release from the endoplasmic reticulum (ER). Increasing intracellular Ca2+ , in turn, provides the excitatory current to the same neuron. Such bidirectional communication can provide long-term potentiation and lead to spontaneous oscillations in the dressed neuron. However, the development of generalized and simplified models for neural–glial ensembles can be useful to study the main types of glion1 response and the resulting dynamical patterns. Such an approach has been widely used to study the dynamical properties of neurons in terms of simplified two-dimensional neuron models 1 Following Ref. Rojtbak (1993) we refer a glial cell as “a glion”. This corresponds to the way “neuron” is used for nerve cells.
0303-2647/$ – see front matter © 2006 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.biosystems.2006.04.012
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
85
(Moris–Lecar, Hindmarch–Rose, FitzHugh–Nagumo). In this paper, the generalized model of a tripartite synapse (Araquo et al., 1999) (i.e., a synaptic connection involving a presynaptic neuron, a postsynaptic neuron, and a glion) is considered. We address the modeling problem using a functional approach that combines the following steps: (i) Temporal and spatial scales should be defined from experimental data; (ii) The important causal chains should be identified from physiological data to be incorporated in the model before the equations are written down; (iii) Simple, well-known and tested functional subunits should be used to model the system. Such a functional model cannot provide the quantitative description of the considered processes because variables and parameters are dimensionless and normalized, and the time scales are fixed as control parameters. However, the model should be able to reproduce the main dynamical patterns of the neural–glial ensemble and allow us to predict its changes with varying control parameters. The paper is organized as follows. In Section 2, we discuss the main pathways of neural–glial interaction and select the facts to be modeled. In Section 3, the mathematical model and the control parameters are described. In Section 4, numerical simulations mimic the glion activation pathways and predict the long-term potentiation effect for the generalized tripartite synapse.
Fig. 1. Main pathways of glion activation (indicated with empty-head arrows) and response (indicated with full-head arrows). Two mechanisms of glion activation are presented: (i) fast mechanism via glion depolarization due to rise of the extracellular potassium concentration (thick solid line) and (ii) slow mechanism associated with the IP3 production (dotted line).
•
• •
2. Neural–glial interaction pathways Let us briefly list the main facts and causal chains related to dynamical activity of the simple neural–glial ensemble to be modeled. The detailed description of biochemical mechanisms can be found, for example, in Refs. Rojtbak (1993), Carmignoto (2000) and Deitmer et al. (1999): • Glial cells do not generate action potentials but exhibit their calcium excitability and self-sustained oscillatory behavior with typical period from 50 ms to 1 min; • Glial cells monitor activity-dependent changes in the chemical environment of the extracellular space shared with neurons. Variation of the potassium concentration is most significant because the intracellular concentration of [K+ ] is much higher than the extracellular one. For example, in ganglia of the medical leech, the extracellular potassium concentration can fluctuate from a normal level of 4 to 10 mM under
• •
conditions of high neuronal activity (Deitmer et al., 1999); A rise of the extracellular [K+ ] instantly causes glion depolarization. Experimentally recorded changes of transmembrane potential are in a good agreement with calculated changes of the potassium equilibrium potential; Glial cells have voltage gated Ca2+ channels that provide an influx of calcium into the cytoplasm in response to glion depolarization; When the synaptic terminal is activated by the presynaptic neuron, some amount of mediator (e.g., glutamate) can leave the synaptic cleft, reach receptors on the glion membrane and, in turn, activate the secondary messenger production (IP3 ). Increasing [IP3 ] triggers the calcium release from the ER and can, thus, evoke Ca2+ oscillations; Increasing [Ca2+ ] in the glion cytoplasm triggers the production of glion mediator (glutamate) and its release into the intercellular space; The glion mediator affects the synaptic terminal providing [Ca2+ ] elevation and depolarization of the postsynaptic neuron. At the same time, the glion mediator causes inhibition of the postsynaptic potentials, i.e., the glion activity to some extent reduces the synaptic strength.
To clarify the structure of the mathematical model that we introduce in the next section, Fig. 1 shows schematically the above listed pathways. Our description of the neural–glial interaction pathways is clearly simplified: several mediators may be involved, increas-
86
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
ing intercellular [K+ ] triggers the removal of potassium by glions, etc. However, in the framework of our approach we take into account four main causal chains. Namely, we consider a fast and a slow mechanisms of glion activation. The fast mechanism is represented by the glion depolarization due to rise of the extracellular potassium concentration. The slow mechanism is associated with the IP3 production initiated by synapse mediator diffusion. We also account for the glial response both to the presynaptic neuron (its synaptic terminal) and to the postsynaptic neuron. 3. Model
ss (v1 − hs ) positive and switches the hyperbolic tangent function to positive values. As a result, z increases to z ≈ 1 with the rate proportional to 1/τs . When v1 is reduced, z inactivates back to zero. Once activated, z provides the postsynaptic current Isyn applied to the second neuron. The factor ks plays the role of conductivity while δGm reflects glion response. Reference level z0 is calculated from the assumption that when the presynaptic neuron is silent then z(t) = z0 . This gives z0 =
2ds . 1 + 2ds + exp(2I1 )
(4)
3.3. Postsynaptic neuron
Although our model is qualitative and dimensionless, we will preserve the widely used notions for the main parameters and variables. In this way, v1 and v2 for the fast variables of the neuron model will represent the transmembrane potentials, Iapp plays the role of the applied current, etc. This will allow us to discuss the model in terms of the processes that we aim to simulate.
The postsynaptic neuron is also described by a FitzHugh–Nagumo model: dv2 v3 = v2 − 2 − w2 , dt 3 dw2 = v2 + I2 − Isyn − Iglion , dt
ε2
(5)
where v1 and w1 are the fast and the slow variables, ε1 = 0.04 is the time separation parameter, I1 = 1.05 is the control parameter. For adjusting the operating regime of the presynaptic neuron, an excitatory current Iapp is added.
with I2 = 1.02 and ε2 = 0.04. Instead of the control parameter Iapp as used in the presynaptic neuron, both the synaptic current Isyn and the glion-induced current Iglion are included in Eq. (5). The glion-induced current is proportional to the glion mediator production Gm with a factor γ, Iglion = γGm . Thus, the parameters δ and γ control the strength of the glion influence on the synapse and the postsynaptic neuron, respectively. In the glion, the processes to be modeled are the slow IP3 production in response to synapse activity, [Ca2+ ] oscillations, and the production Gm of glion mediator as controlled by the calcium concentration.
3.2. Synapse
3.4. Calcium dynamics of glion
In the framework of our approach, the essential properties of synaptic coupling between two neurons are the delayed response of the postsynaptic neuron activity, and the threshold activation (sigmoid nonlinearity). Following Kopell et al. (2000), we describe the synaptic coupling with the first-order differential equation:
Since we use the FitzHugh–Nagumo equations for the two neurons, it is reasonably to adopt a simple and qualitative model for the intra-glion calcium dynamics. Hence, we apply the two-dimensional two-pool model (Keener and Sneyd, 1998) with additional terms describing the external forcing:
3.1. Presynaptic neuron This neuron is described by the well-known FitzHugh–Nagumo model (FitzHugh, 1961): ε1
τs
dv1 v3 = v1 − 1 − w1 , dt 3
dw1 = v1 + I1 − Iapp , dt
dz z = (1 + tanh(ss (v1 − hs )))(1 − z) − , dt ds
Isyn = (ks − δGm )(z − z0 ).
(1)
(2) (3)
Here, z is a synaptic activation variable and τs describes the time delay. The parameters hs , ss , and ds are responsible for activation and relaxation of z. When v1 < hs , the synapse is inactive and z ≈ 0. Increasing v1 makes
τc
dc = −c − c4 f (c, ce ) + (r + αw2 + βSm ), dt
εc τc
dce = f (c, ce ). dt
(6) (7)
Here, c denotes the Ca2+ concentration within the glion, ce is the calcium concentration in the internal store (endoplasmic reticulum ER), τc defines the characteristic
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
time for Ca2+ oscillations together with time separation parameter εc . The value of the parameter r = 0.31 controls the initial state of the calcium oscillator without external influence (α = 0 and β = 0). The term (r + αw2 + βSm ) represents the calcium influx from the extracellular space that is sensitive to the synapse mediator Sm production (with the factor β) and to the glion depolarization by increasing extracellular potassium (via αw2 ). The Ca2+ exchange between the cytoplasm and ER is described with a nonlinear function (Keener and Sneyd, 1998): 2 c4 c2 ce −c3 ce . f (c, ce ) = c1 − 1 + c2 1 + ce2 c24 + c4 (8) 3.5. Mediator production Finally, we take into account the secondary mediator Sm (IP3 ) and glion mediator Gm productions. Both processes are characterized by the considerable delay. There are threshold value for the Sm production that is triggered by synaptic activity (z variable) and for the Gm release that is governed by the rise of the calcium concentration in the glion. At qualitative level, we assume that the dynamics of Sm and Gm are governed by similar equations: τSm
dSm Sm =(1 + tanh(sSm (z − hSm ))) × (1 − Sm ) − , dt dSm
87
the processes under modeling. In this way, the reference value for all time scales is equal to 1.0. It is given by the second equation for both neurons. Thus, one can adjust the synaptic delay by varying τs ∈ [0.5; 10.0]. According to our assumptions the following relation between the other time scales has to be maintained: τSm τGm τs .
(11)
Threshold parameters hs , hSm , and hGm are selected to distinguish between inactivated and activated states of the variables v1 , z, and c, respectively. Parameters ds , dSm , and dGm at fixed at some typical values estimated from Ref. Kopell et al. (2000). Below we use α, β, γ, and δ as the main control parameters involved in the regulation of the glion activation and response. Iapp is used as a variable external current that triggers spontaneous activity of the presynaptic neuron and, thus, allows us to consider the transient dynamics of the neural–glial ensemble. Let us consider the functional structure of the model displayed in Fig. 2. It is clear that only two types of subunits are involved in the description: oscillatory/excitable subunits and threshold-activated subunits. Note that in spite of the neuron Eqs. (1) and (5), and the calcium model (6) are different in form, they belong to the same type of two-dimensional oscillators and, thus, (6) can be converted into the form of the FitzHugh– Nagumo model with the appropriate substitution of variables (Keener and Sneyd, 1998). This structure
(9) τGm
dGm = (1 + tanh(sGm (c − hGm ))) × (1 − Gm ) dt Gm , (10) − dGm
where the specific values of the time scales, threshold values, and steepnesses of the sigmoid functions are determined by an appropriate choice of the control parameters τSm , τGm , sSm , sGm , hSm , hGm , dSm , and dGm . A set of control parameters is shown in the following table: c1 = 0.13, c4 = 2.0/εc , τc = 8.0 ss = 1.0, hs = 0.0, ds = 3.0,
c2 = 0.9, τSm = 100.0, sSm = 100.0, hSm = 0.45, dSm = 3.0,
c3 = 0.004, εc = 0.04, τGm = 5.0τs sGm = 100.0 hGm = 0.5 dGm = 3.0
The first six parameters (two top rows) describing the calcium dynamics are taken from Ref. Keener and Sneyd (1998) while remaining parameters were selected to comply qualitatively with the known relations between
Fig. 2. Functional diagram of the subunits with corresponding variables of the model. Both the two neurons and the glion are oscillatory/ excitable units while synapse and glion mediator productions are modeled using threshold-activation units. The control parameters α and β are related to the fast (thick solid line) and slow (dotted line) mechanisms of activation, respectively, while parameters δ and γ specify glion response.
88
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
allows us (i) to adjust the model to the specific problem by an appropriate selection of time scales/threshold values and (ii) to develop the simplest possible paradigm for a neural–glial ensemble like integrate-and-fire model being the simplest neuron paradigm. 4. Dynamical properties of neuron–glion interaction 4.1. Dynamics induced via the fast activation pathway To test the model we block the glion response by setting γ = 0 and δ = 0. With this parameter choice, the glion can be activated by the neurons, but there is no feedback to the neural dynamics. If β = 0, the slow activation pathway is disconnected. To consider the transient dynamics of the intra-glion calcium concentration c, the presynaptic neuron is stimulated by Iapp during some time interval. With synapse parameters at ks = 0.2 and τs = 10.0, the postsynaptic neuron also generates spike train with a time delay defined by τs . The second variable w2 representing the potassium current of the postsynaptic neuron, affects the calcium oscillator of the glion with a factor α (Eq. (6)). Fig. 3 shows the temporal dynamics of the transmembrane potential variable v2 together with representative time courses for c at two different values of α. At α = 0.01, the activity of the postsynaptic neuron evokes subthreshold damped oscillations (Fig. 3(a)). Note that subthreshold response is initiated at the beginning of the spike train (increasing c) and at the end of the spike train (decreasing c) but not during the firing. With α = 0.02 (Fig. 3(b)), the forcing is strong enough to evoke a single spike of calcium oscillator at the beginning of neuron spike train. Thus, the fast activation pathway leads to excitable dynamics and seems to be sensitive to sharp changes of activity in the postsynaptic neuron. The slow calcium oscillator just averages the fast neuron firings and reacts on changes of the mean value of w2 . Note that there is the rapid activation of subthreshold oscillations, but relaxation to the unperturbed state is relatively slow. Like resonator-type systems (Izhikevich, 2001), whose response strongly depends on the forcing frequency, the calcium oscillator can produce more complicated dynamical patterns in response to a complex bursting activity of the neuron.
Fig. 3. The fast activation pathway of glion leads to a calcium resonator-type excitable dynamics. The top panel represents the dynamics of the postsynaptic neuron. The two bottom panels show temporal variations of calcium concentration in the glion: (a) subthreshold response at α = 0.01 and (b) single-spike response at α = 0.02. Note the different scale for c variable in (a) and (b).
via secondary mediator production with varying β (Eqs. (2) and (9)). Because τSm τc , the secondary mediator production is much slower than the calcium oscillations. Thus, both the factor β and the duration of neuron firing play a role in the activation of the glion. Let us first consider the influence of β at the fixed firing interval of the presynaptic neuron. The time course of v1 and the response of the calcium oscillator with increasing β are shown in Fig. 4. At β = 0.003 (Fig. 4(a)), there are subthreshold calcium oscillations in the glion. Their amplitude slowly increases during the period of presynaptic neuron activity. When the presynaptic neuron stops firing, the amplitude of the subthreshold oscillations still increases to finally gradually decrease.
4.2. Patterns generated by slow activation pathway Let us now block the fast activation pathway by setting α = 0 but allow a synapse to activate the glion
Fig. 4. Slow pathway of glion activation. The presynaptic neuron activity (top panel) evokes different glion responses (a–d) at different rates of secondary messenger production β.
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
Fig. 5. Depending on the duration T of presynaptic neuron activity, the glion shows different responses of its calcium dynamics: (a) At T = 1000, time interval is too short to evoke full-scale oscillations in the glion, and only subthreshold oscillations are observed. At T = 2000 (b) and T = 3000 (c), extended activity of the presynaptic neuron evokes a self-sustained response. Note the different scale for variable c in (a–c).
The calcium oscillator does not produce any spikes, and the glion remains inactivated. At higher β = 0.0033 (Fig. 4(b)), the amplitude increases faster, and a single full-scale spike is produced at the end of the firing course of the presynaptic neuron. This resembles the excitatory dynamics shown in Fig. 3. However, the further increase of β (Fig. 4(c) and (d)) reveals another mechanism. Spikes occur closer and closer to the beginning of firing course implying faster activation of the glion. Let us now set a smaller value of β = 0.002 but increase considerably the duration of the presynaptic neuron firing T (Fig. 5). At T = 1000 (a), a subthreshold response can be observed. This is consistent with the previous figure. At T = 2000 (b), the calcium oscillator starts spiking at t ≈ 1750 but stops immediately after the presynaptic neuron stops its firing. With further increase of T to 2000 (c), calcium spiking is maintained during the neuron firing but it is terminated as soon as the neuron becomes silent. The interspike intervals seem to be approximately the same, and no pronounced frequency variation is observed for small β. Note, the glion inactivation time at the considered values of β is quite short: as soon as the neuron stops firing, the calcium oscillator also stops producing spikes. What will happen at large β = 0.5 when external forcing of the calcium oscillator is strong? It is clearly seen in Fig. 6 that the response of the calcium oscillator is more complicated. At the beginning of the neuron firing period, there is a short burst of activity, then the calcium
89
Fig. 6. At large rate β of secondary messenger Sm production, the calcium oscillations in the glion can be blocked. The bottom panel shows that this is caused by the high level of Sm .
oscillations become blocked and the calcium concentration c reaches some constant level higher than in the inactivated state. When the stimulation is terminated, calcium oscillations with increasing interspike intervals are evoked. Finally, the calcium oscillator returns to its resting state. The last panel shows the rapid increase of Sm , a plateau with constant value, and finally a gradual fall. Obviously, only intermediate values of Sm evoke calcium oscillations, and not the maximum values. The observed behavior can again be explained in terms of the dynamics of the calcium oscillator. Fig. 7 shows the bifurcation diagram for calcium subsystem (Eq. (6) with α = 0, β = 0). There are two Andronov–Hopf bifurcations at r1 and r2 . Between these points the system exhibits oscillations whose amplitude decreases while the frequency increases with increasing r. Thus, the transient dynamics observed in Fig. 6 can be explained by a shift of operating point: increasing term βSm corresponds to large values of r and the calcium subsystem can be driven out the oscillatory regime.
Fig. 7. One-parameter bifurcational diagram for calcium oscillator (Eq. (6)) indicates the limited range of oscillatory behavior. During neuronal activity, the operating point can be shifted along r from one equilibrium state to another.
90
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
that to make clear pictures of these patterns, we have used τSm = 100. This value satisfies Eq. (11), but realistic values could be larger. 5. Conclusion
Fig. 8. Long-term potentiation. Top panels: temporal dynamics of the presynaptic neuron activity v1 , glion calcium concentration c, and glion-produced mediator Gm . Bottom panels: evoked activity of the postsynaptic neuron at different γ.
4.3. Long-term potentiation From the above results, we have learned how the glion Ca2+ system can be activated by neurons. Let us now consider the closed loop of neural–glial activation and response to see how the glion can control the activity of the postsynaptic neuron. The control parameters are fixed at α = 0, β = 0.05, ks = 0.2, τs = 10.0, and δ = 0 when the fast activation pathway as well as the glion-induced inactivation of the synapse are blocked. The top panels in Fig. 8 show the glion activation c and corresponding production of mediator Gm that are almost twice as long as the duration of the presynaptic neuron firing. The bottom three panels in Fig. 8 represent the response of the postsynaptic neuron at different values of γ. At γ = 0.059, there is a single additional burst when the presynaptic neuron is already in the resting state. This burst is observed at the moment of peak in the Gm production. However, the next Gm peaks of almost the same amplitude do not produce neuronal firing. Only a subthreshold increase of v2 can be observed. At γ = 0.07, the glial response produces a sequence of bursts with decreasing duration. Thus, the glion evokes complex firing patterns that are different from the excitatory signal produced by the presynaptic neuron. Finally, at γ = 0.5, a continuous firing of the postsynaptic neuron is maintained during all the time interval of Gm production. We did not observe any bursting at a tail of the response. We conclude that varying γ produces different patterns of potentiation of the postsynaptic neuron. Note
Based on the simplified description of neural–glial interaction pathways in a tripartite synapse we proposed the generalized model (1)–(10) that preserves essential features of such a functional unit: fast and slow activation pathways as well as a dual response of the glion to synaptic and neuronal activities. The model dynamics clearly resembles experimental observations. Namely, there are two types of glion activation (excitable dynamics and self-sustained oscillations) as described by Rojtbak (1993). We did not model the particular receptors and their activation but the response of Ca2+ oscillator to the external signal seems to be generic and, thus, our model reproduces it at fast and slow time scales. Our model clearly predicts the effect of long-term potentiation that is believed to be the main feature of glial cells. It is interesting that there are different types of post-excitatory behavior, not just continuous firing. Such model prediction needs experimental verification on specific objects. Since our model is developed to simulate dynamics of neural–glial ensembles, its extension to multi-synapse and multi-neuron structures is of great interest. This can be done by modifying of Eq. (6) to the following form: ⎞ ⎛ N M dc τc = −c−c4 f (c, ce )+⎝r+ αi w i + βj Smj ⎠, dt i=1
εc τc
dce = f (c, ce ), dt
j=1
(12)
where i = 1, . . . , N enumerates the synaptic terminals, and j = 1, . . . , M counts the neurons that are closely located to the glion (generally, N = M). In Eq. (12), c can represent not only single-cell calcium dynamics but also its variation in small networks of glions coupled via gap junctions. Acknowledgments This work was supported by the European Union through the Network of Excellence BioSim (Contract No. LSHB-CT-2004-005137) and RFBR grant 04-0216769. The authors acknowledge fruitful discussions with N. Brazhe, A. Brazhe, and G. Maksimov (Moscow State University).
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
References Araquo, A., Parpura, V., Zanzgiri, R.P., Haydon, P.H., 1999. Tripartite synapses: glia, the unacknowledged. Trends Neurosci. 22, 208. Bonvento, G., Giaume, C., Lorenceau, J., 2002. Neuron–glia interactions: from physiology to behavior. J. Physiol. 96, 167– 168. Carmignoto, G., 2000. Reciprocal communication systems between astrocytes and neurons. Prog. Neurobiol. 62, 561–581. Deitmer, J.W., Rose, C.R., Munch, T., Schmidt, J., Nett, W., Schneider, N.-P., Lohr, C., 1999. Leech glial cell: functional role in a simple nervous system. Glia 28, 175–182. Fields, R.D., Stevens-Graham, B., 2002. New insights into neuron–glia communication. Science 298, 556–562.
91
FitzHugh, R.A., 1961. Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–446. Izhikevich, E.M., 2001. Resonate-and-fire neurons. Neural Netw. J. 14, 883–894. Keener, J., Sneyd, J., 1998. Mathematical Physiology. Springer, New York. Kopell, N., Ermentrout, G.B., Whittington, M.A., Traub, R.D., 2000. Gamma rhythms and beta rhythms have different synchronization properties. Proc. Natl. Acad. Sci. U.S.A. 97, 1867–1872. Nadkarni, S., Jung, P., 2003. Spontaneous oscillations of dressed neurons: a new mechanism for epilepsy?. Phys. Rev. Lett. 91, 268101. Nadkarni, S., Jung, P, 2004. Dressed neurons: modeling neural–glia interactions. Phys. Biol. 1, 35–41. Rojtbak, A.I., 1993. Glia. Nauka, Sankt-Petersburg (in Russian).
Functional modeling of neural–glial interaction D.E. Postnova , L.S. Ryazanovaa , O.V. Sosnovtsevab,∗ a
Physics Department, Saratov State University, Astrakhanskaya Street 83, Saratov 410026, Russia Department of Physics, The Technical University of Denmark, 2800 Kongens Lyngby, Denmark
b
Received 24 December 2005; accepted 15 April 2006
Abstract We propose a generalized mathematical model for a small neural–glial ensemble. The model incorporates subunits of the tripartite synapse that includes a presynaptic neuron, the synaptic terminal itself, a postsynaptic neuron, and a glial cell. The glial cell is assumed to be activated via two different pathways: (i) the fast increase of intercellular [K+ ] produced by the spiking activity of the postsynaptic neuron, and (ii) the slow production of a mediator triggered by the synaptic activity. Our model predicts the long-term potentiation of the postsynaptic neuron as well as various [Ca2+ ] transients in response to the activation of different pathways. © 2006 Elsevier Ireland Ltd. All rights reserved. Keywords: Glion; Neural–glial interaction; Functional modeling
1. Introduction During the last decade, application of a variety of new tools has allowed us to uncover the diverse dynamic roles of glial cells. Glial cells are active partners to neurons in processing information and synaptic integration (Field and Steven-Graham, 2002; Bonvento et al., 2002). Thus, studying the interaction between the neurons and glial cells has become an important problem in neurophysiology and cell biophysics. The modeling of neural–glial interaction may help to understand the transition between normal and different pathological states. An adequate mathematical model of neural–glial interaction should account for all relevant ionic currents, processes of production, propagation of neuromediators, etc. Such a model inevitably will be high-dimensional with a large number of control parameters. While the
∗
Corresponding author. Tel.: +45 4525 3106; fax: +45 4593 1669. E-mail addresses: [email protected] (D.E. Postnov), [email protected] (L.S. Ryazanova), [email protected] (O.V. Sosnovtseva). URL: http://www.fys.dtu.dk/∼olga.
main pathways of neural–glial interactions seem to be quite common, the details can vary from species to species. Hence, there will be a number of speciesoriented mathematical models. Nadkarni and Jung (2003, 2004) first introduced the concept of a “dressed neuron”, representing a set of equations that describe the activation of an astrocyte by a neuron via the inositol 1,4,5-triphosphate (IP3 ) production that triggers the calcium release from the endoplasmic reticulum (ER). Increasing intracellular Ca2+ , in turn, provides the excitatory current to the same neuron. Such bidirectional communication can provide long-term potentiation and lead to spontaneous oscillations in the dressed neuron. However, the development of generalized and simplified models for neural–glial ensembles can be useful to study the main types of glion1 response and the resulting dynamical patterns. Such an approach has been widely used to study the dynamical properties of neurons in terms of simplified two-dimensional neuron models 1 Following Ref. Rojtbak (1993) we refer a glial cell as “a glion”. This corresponds to the way “neuron” is used for nerve cells.
0303-2647/$ – see front matter © 2006 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.biosystems.2006.04.012
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
85
(Moris–Lecar, Hindmarch–Rose, FitzHugh–Nagumo). In this paper, the generalized model of a tripartite synapse (Araquo et al., 1999) (i.e., a synaptic connection involving a presynaptic neuron, a postsynaptic neuron, and a glion) is considered. We address the modeling problem using a functional approach that combines the following steps: (i) Temporal and spatial scales should be defined from experimental data; (ii) The important causal chains should be identified from physiological data to be incorporated in the model before the equations are written down; (iii) Simple, well-known and tested functional subunits should be used to model the system. Such a functional model cannot provide the quantitative description of the considered processes because variables and parameters are dimensionless and normalized, and the time scales are fixed as control parameters. However, the model should be able to reproduce the main dynamical patterns of the neural–glial ensemble and allow us to predict its changes with varying control parameters. The paper is organized as follows. In Section 2, we discuss the main pathways of neural–glial interaction and select the facts to be modeled. In Section 3, the mathematical model and the control parameters are described. In Section 4, numerical simulations mimic the glion activation pathways and predict the long-term potentiation effect for the generalized tripartite synapse.
Fig. 1. Main pathways of glion activation (indicated with empty-head arrows) and response (indicated with full-head arrows). Two mechanisms of glion activation are presented: (i) fast mechanism via glion depolarization due to rise of the extracellular potassium concentration (thick solid line) and (ii) slow mechanism associated with the IP3 production (dotted line).
•
• •
2. Neural–glial interaction pathways Let us briefly list the main facts and causal chains related to dynamical activity of the simple neural–glial ensemble to be modeled. The detailed description of biochemical mechanisms can be found, for example, in Refs. Rojtbak (1993), Carmignoto (2000) and Deitmer et al. (1999): • Glial cells do not generate action potentials but exhibit their calcium excitability and self-sustained oscillatory behavior with typical period from 50 ms to 1 min; • Glial cells monitor activity-dependent changes in the chemical environment of the extracellular space shared with neurons. Variation of the potassium concentration is most significant because the intracellular concentration of [K+ ] is much higher than the extracellular one. For example, in ganglia of the medical leech, the extracellular potassium concentration can fluctuate from a normal level of 4 to 10 mM under
• •
conditions of high neuronal activity (Deitmer et al., 1999); A rise of the extracellular [K+ ] instantly causes glion depolarization. Experimentally recorded changes of transmembrane potential are in a good agreement with calculated changes of the potassium equilibrium potential; Glial cells have voltage gated Ca2+ channels that provide an influx of calcium into the cytoplasm in response to glion depolarization; When the synaptic terminal is activated by the presynaptic neuron, some amount of mediator (e.g., glutamate) can leave the synaptic cleft, reach receptors on the glion membrane and, in turn, activate the secondary messenger production (IP3 ). Increasing [IP3 ] triggers the calcium release from the ER and can, thus, evoke Ca2+ oscillations; Increasing [Ca2+ ] in the glion cytoplasm triggers the production of glion mediator (glutamate) and its release into the intercellular space; The glion mediator affects the synaptic terminal providing [Ca2+ ] elevation and depolarization of the postsynaptic neuron. At the same time, the glion mediator causes inhibition of the postsynaptic potentials, i.e., the glion activity to some extent reduces the synaptic strength.
To clarify the structure of the mathematical model that we introduce in the next section, Fig. 1 shows schematically the above listed pathways. Our description of the neural–glial interaction pathways is clearly simplified: several mediators may be involved, increas-
86
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
ing intercellular [K+ ] triggers the removal of potassium by glions, etc. However, in the framework of our approach we take into account four main causal chains. Namely, we consider a fast and a slow mechanisms of glion activation. The fast mechanism is represented by the glion depolarization due to rise of the extracellular potassium concentration. The slow mechanism is associated with the IP3 production initiated by synapse mediator diffusion. We also account for the glial response both to the presynaptic neuron (its synaptic terminal) and to the postsynaptic neuron. 3. Model
ss (v1 − hs ) positive and switches the hyperbolic tangent function to positive values. As a result, z increases to z ≈ 1 with the rate proportional to 1/τs . When v1 is reduced, z inactivates back to zero. Once activated, z provides the postsynaptic current Isyn applied to the second neuron. The factor ks plays the role of conductivity while δGm reflects glion response. Reference level z0 is calculated from the assumption that when the presynaptic neuron is silent then z(t) = z0 . This gives z0 =
2ds . 1 + 2ds + exp(2I1 )
(4)
3.3. Postsynaptic neuron
Although our model is qualitative and dimensionless, we will preserve the widely used notions for the main parameters and variables. In this way, v1 and v2 for the fast variables of the neuron model will represent the transmembrane potentials, Iapp plays the role of the applied current, etc. This will allow us to discuss the model in terms of the processes that we aim to simulate.
The postsynaptic neuron is also described by a FitzHugh–Nagumo model: dv2 v3 = v2 − 2 − w2 , dt 3 dw2 = v2 + I2 − Isyn − Iglion , dt
ε2
(5)
where v1 and w1 are the fast and the slow variables, ε1 = 0.04 is the time separation parameter, I1 = 1.05 is the control parameter. For adjusting the operating regime of the presynaptic neuron, an excitatory current Iapp is added.
with I2 = 1.02 and ε2 = 0.04. Instead of the control parameter Iapp as used in the presynaptic neuron, both the synaptic current Isyn and the glion-induced current Iglion are included in Eq. (5). The glion-induced current is proportional to the glion mediator production Gm with a factor γ, Iglion = γGm . Thus, the parameters δ and γ control the strength of the glion influence on the synapse and the postsynaptic neuron, respectively. In the glion, the processes to be modeled are the slow IP3 production in response to synapse activity, [Ca2+ ] oscillations, and the production Gm of glion mediator as controlled by the calcium concentration.
3.2. Synapse
3.4. Calcium dynamics of glion
In the framework of our approach, the essential properties of synaptic coupling between two neurons are the delayed response of the postsynaptic neuron activity, and the threshold activation (sigmoid nonlinearity). Following Kopell et al. (2000), we describe the synaptic coupling with the first-order differential equation:
Since we use the FitzHugh–Nagumo equations for the two neurons, it is reasonably to adopt a simple and qualitative model for the intra-glion calcium dynamics. Hence, we apply the two-dimensional two-pool model (Keener and Sneyd, 1998) with additional terms describing the external forcing:
3.1. Presynaptic neuron This neuron is described by the well-known FitzHugh–Nagumo model (FitzHugh, 1961): ε1
τs
dv1 v3 = v1 − 1 − w1 , dt 3
dw1 = v1 + I1 − Iapp , dt
dz z = (1 + tanh(ss (v1 − hs )))(1 − z) − , dt ds
Isyn = (ks − δGm )(z − z0 ).
(1)
(2) (3)
Here, z is a synaptic activation variable and τs describes the time delay. The parameters hs , ss , and ds are responsible for activation and relaxation of z. When v1 < hs , the synapse is inactive and z ≈ 0. Increasing v1 makes
τc
dc = −c − c4 f (c, ce ) + (r + αw2 + βSm ), dt
εc τc
dce = f (c, ce ). dt
(6) (7)
Here, c denotes the Ca2+ concentration within the glion, ce is the calcium concentration in the internal store (endoplasmic reticulum ER), τc defines the characteristic
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
time for Ca2+ oscillations together with time separation parameter εc . The value of the parameter r = 0.31 controls the initial state of the calcium oscillator without external influence (α = 0 and β = 0). The term (r + αw2 + βSm ) represents the calcium influx from the extracellular space that is sensitive to the synapse mediator Sm production (with the factor β) and to the glion depolarization by increasing extracellular potassium (via αw2 ). The Ca2+ exchange between the cytoplasm and ER is described with a nonlinear function (Keener and Sneyd, 1998): 2 c4 c2 ce −c3 ce . f (c, ce ) = c1 − 1 + c2 1 + ce2 c24 + c4 (8) 3.5. Mediator production Finally, we take into account the secondary mediator Sm (IP3 ) and glion mediator Gm productions. Both processes are characterized by the considerable delay. There are threshold value for the Sm production that is triggered by synaptic activity (z variable) and for the Gm release that is governed by the rise of the calcium concentration in the glion. At qualitative level, we assume that the dynamics of Sm and Gm are governed by similar equations: τSm
dSm Sm =(1 + tanh(sSm (z − hSm ))) × (1 − Sm ) − , dt dSm
87
the processes under modeling. In this way, the reference value for all time scales is equal to 1.0. It is given by the second equation for both neurons. Thus, one can adjust the synaptic delay by varying τs ∈ [0.5; 10.0]. According to our assumptions the following relation between the other time scales has to be maintained: τSm τGm τs .
(11)
Threshold parameters hs , hSm , and hGm are selected to distinguish between inactivated and activated states of the variables v1 , z, and c, respectively. Parameters ds , dSm , and dGm at fixed at some typical values estimated from Ref. Kopell et al. (2000). Below we use α, β, γ, and δ as the main control parameters involved in the regulation of the glion activation and response. Iapp is used as a variable external current that triggers spontaneous activity of the presynaptic neuron and, thus, allows us to consider the transient dynamics of the neural–glial ensemble. Let us consider the functional structure of the model displayed in Fig. 2. It is clear that only two types of subunits are involved in the description: oscillatory/excitable subunits and threshold-activated subunits. Note that in spite of the neuron Eqs. (1) and (5), and the calcium model (6) are different in form, they belong to the same type of two-dimensional oscillators and, thus, (6) can be converted into the form of the FitzHugh– Nagumo model with the appropriate substitution of variables (Keener and Sneyd, 1998). This structure
(9) τGm
dGm = (1 + tanh(sGm (c − hGm ))) × (1 − Gm ) dt Gm , (10) − dGm
where the specific values of the time scales, threshold values, and steepnesses of the sigmoid functions are determined by an appropriate choice of the control parameters τSm , τGm , sSm , sGm , hSm , hGm , dSm , and dGm . A set of control parameters is shown in the following table: c1 = 0.13, c4 = 2.0/εc , τc = 8.0 ss = 1.0, hs = 0.0, ds = 3.0,
c2 = 0.9, τSm = 100.0, sSm = 100.0, hSm = 0.45, dSm = 3.0,
c3 = 0.004, εc = 0.04, τGm = 5.0τs sGm = 100.0 hGm = 0.5 dGm = 3.0
The first six parameters (two top rows) describing the calcium dynamics are taken from Ref. Keener and Sneyd (1998) while remaining parameters were selected to comply qualitatively with the known relations between
Fig. 2. Functional diagram of the subunits with corresponding variables of the model. Both the two neurons and the glion are oscillatory/ excitable units while synapse and glion mediator productions are modeled using threshold-activation units. The control parameters α and β are related to the fast (thick solid line) and slow (dotted line) mechanisms of activation, respectively, while parameters δ and γ specify glion response.
88
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
allows us (i) to adjust the model to the specific problem by an appropriate selection of time scales/threshold values and (ii) to develop the simplest possible paradigm for a neural–glial ensemble like integrate-and-fire model being the simplest neuron paradigm. 4. Dynamical properties of neuron–glion interaction 4.1. Dynamics induced via the fast activation pathway To test the model we block the glion response by setting γ = 0 and δ = 0. With this parameter choice, the glion can be activated by the neurons, but there is no feedback to the neural dynamics. If β = 0, the slow activation pathway is disconnected. To consider the transient dynamics of the intra-glion calcium concentration c, the presynaptic neuron is stimulated by Iapp during some time interval. With synapse parameters at ks = 0.2 and τs = 10.0, the postsynaptic neuron also generates spike train with a time delay defined by τs . The second variable w2 representing the potassium current of the postsynaptic neuron, affects the calcium oscillator of the glion with a factor α (Eq. (6)). Fig. 3 shows the temporal dynamics of the transmembrane potential variable v2 together with representative time courses for c at two different values of α. At α = 0.01, the activity of the postsynaptic neuron evokes subthreshold damped oscillations (Fig. 3(a)). Note that subthreshold response is initiated at the beginning of the spike train (increasing c) and at the end of the spike train (decreasing c) but not during the firing. With α = 0.02 (Fig. 3(b)), the forcing is strong enough to evoke a single spike of calcium oscillator at the beginning of neuron spike train. Thus, the fast activation pathway leads to excitable dynamics and seems to be sensitive to sharp changes of activity in the postsynaptic neuron. The slow calcium oscillator just averages the fast neuron firings and reacts on changes of the mean value of w2 . Note that there is the rapid activation of subthreshold oscillations, but relaxation to the unperturbed state is relatively slow. Like resonator-type systems (Izhikevich, 2001), whose response strongly depends on the forcing frequency, the calcium oscillator can produce more complicated dynamical patterns in response to a complex bursting activity of the neuron.
Fig. 3. The fast activation pathway of glion leads to a calcium resonator-type excitable dynamics. The top panel represents the dynamics of the postsynaptic neuron. The two bottom panels show temporal variations of calcium concentration in the glion: (a) subthreshold response at α = 0.01 and (b) single-spike response at α = 0.02. Note the different scale for c variable in (a) and (b).
via secondary mediator production with varying β (Eqs. (2) and (9)). Because τSm τc , the secondary mediator production is much slower than the calcium oscillations. Thus, both the factor β and the duration of neuron firing play a role in the activation of the glion. Let us first consider the influence of β at the fixed firing interval of the presynaptic neuron. The time course of v1 and the response of the calcium oscillator with increasing β are shown in Fig. 4. At β = 0.003 (Fig. 4(a)), there are subthreshold calcium oscillations in the glion. Their amplitude slowly increases during the period of presynaptic neuron activity. When the presynaptic neuron stops firing, the amplitude of the subthreshold oscillations still increases to finally gradually decrease.
4.2. Patterns generated by slow activation pathway Let us now block the fast activation pathway by setting α = 0 but allow a synapse to activate the glion
Fig. 4. Slow pathway of glion activation. The presynaptic neuron activity (top panel) evokes different glion responses (a–d) at different rates of secondary messenger production β.
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
Fig. 5. Depending on the duration T of presynaptic neuron activity, the glion shows different responses of its calcium dynamics: (a) At T = 1000, time interval is too short to evoke full-scale oscillations in the glion, and only subthreshold oscillations are observed. At T = 2000 (b) and T = 3000 (c), extended activity of the presynaptic neuron evokes a self-sustained response. Note the different scale for variable c in (a–c).
The calcium oscillator does not produce any spikes, and the glion remains inactivated. At higher β = 0.0033 (Fig. 4(b)), the amplitude increases faster, and a single full-scale spike is produced at the end of the firing course of the presynaptic neuron. This resembles the excitatory dynamics shown in Fig. 3. However, the further increase of β (Fig. 4(c) and (d)) reveals another mechanism. Spikes occur closer and closer to the beginning of firing course implying faster activation of the glion. Let us now set a smaller value of β = 0.002 but increase considerably the duration of the presynaptic neuron firing T (Fig. 5). At T = 1000 (a), a subthreshold response can be observed. This is consistent with the previous figure. At T = 2000 (b), the calcium oscillator starts spiking at t ≈ 1750 but stops immediately after the presynaptic neuron stops its firing. With further increase of T to 2000 (c), calcium spiking is maintained during the neuron firing but it is terminated as soon as the neuron becomes silent. The interspike intervals seem to be approximately the same, and no pronounced frequency variation is observed for small β. Note, the glion inactivation time at the considered values of β is quite short: as soon as the neuron stops firing, the calcium oscillator also stops producing spikes. What will happen at large β = 0.5 when external forcing of the calcium oscillator is strong? It is clearly seen in Fig. 6 that the response of the calcium oscillator is more complicated. At the beginning of the neuron firing period, there is a short burst of activity, then the calcium
89
Fig. 6. At large rate β of secondary messenger Sm production, the calcium oscillations in the glion can be blocked. The bottom panel shows that this is caused by the high level of Sm .
oscillations become blocked and the calcium concentration c reaches some constant level higher than in the inactivated state. When the stimulation is terminated, calcium oscillations with increasing interspike intervals are evoked. Finally, the calcium oscillator returns to its resting state. The last panel shows the rapid increase of Sm , a plateau with constant value, and finally a gradual fall. Obviously, only intermediate values of Sm evoke calcium oscillations, and not the maximum values. The observed behavior can again be explained in terms of the dynamics of the calcium oscillator. Fig. 7 shows the bifurcation diagram for calcium subsystem (Eq. (6) with α = 0, β = 0). There are two Andronov–Hopf bifurcations at r1 and r2 . Between these points the system exhibits oscillations whose amplitude decreases while the frequency increases with increasing r. Thus, the transient dynamics observed in Fig. 6 can be explained by a shift of operating point: increasing term βSm corresponds to large values of r and the calcium subsystem can be driven out the oscillatory regime.
Fig. 7. One-parameter bifurcational diagram for calcium oscillator (Eq. (6)) indicates the limited range of oscillatory behavior. During neuronal activity, the operating point can be shifted along r from one equilibrium state to another.
90
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
that to make clear pictures of these patterns, we have used τSm = 100. This value satisfies Eq. (11), but realistic values could be larger. 5. Conclusion
Fig. 8. Long-term potentiation. Top panels: temporal dynamics of the presynaptic neuron activity v1 , glion calcium concentration c, and glion-produced mediator Gm . Bottom panels: evoked activity of the postsynaptic neuron at different γ.
4.3. Long-term potentiation From the above results, we have learned how the glion Ca2+ system can be activated by neurons. Let us now consider the closed loop of neural–glial activation and response to see how the glion can control the activity of the postsynaptic neuron. The control parameters are fixed at α = 0, β = 0.05, ks = 0.2, τs = 10.0, and δ = 0 when the fast activation pathway as well as the glion-induced inactivation of the synapse are blocked. The top panels in Fig. 8 show the glion activation c and corresponding production of mediator Gm that are almost twice as long as the duration of the presynaptic neuron firing. The bottom three panels in Fig. 8 represent the response of the postsynaptic neuron at different values of γ. At γ = 0.059, there is a single additional burst when the presynaptic neuron is already in the resting state. This burst is observed at the moment of peak in the Gm production. However, the next Gm peaks of almost the same amplitude do not produce neuronal firing. Only a subthreshold increase of v2 can be observed. At γ = 0.07, the glial response produces a sequence of bursts with decreasing duration. Thus, the glion evokes complex firing patterns that are different from the excitatory signal produced by the presynaptic neuron. Finally, at γ = 0.5, a continuous firing of the postsynaptic neuron is maintained during all the time interval of Gm production. We did not observe any bursting at a tail of the response. We conclude that varying γ produces different patterns of potentiation of the postsynaptic neuron. Note
Based on the simplified description of neural–glial interaction pathways in a tripartite synapse we proposed the generalized model (1)–(10) that preserves essential features of such a functional unit: fast and slow activation pathways as well as a dual response of the glion to synaptic and neuronal activities. The model dynamics clearly resembles experimental observations. Namely, there are two types of glion activation (excitable dynamics and self-sustained oscillations) as described by Rojtbak (1993). We did not model the particular receptors and their activation but the response of Ca2+ oscillator to the external signal seems to be generic and, thus, our model reproduces it at fast and slow time scales. Our model clearly predicts the effect of long-term potentiation that is believed to be the main feature of glial cells. It is interesting that there are different types of post-excitatory behavior, not just continuous firing. Such model prediction needs experimental verification on specific objects. Since our model is developed to simulate dynamics of neural–glial ensembles, its extension to multi-synapse and multi-neuron structures is of great interest. This can be done by modifying of Eq. (6) to the following form: ⎞ ⎛ N M dc τc = −c−c4 f (c, ce )+⎝r+ αi w i + βj Smj ⎠, dt i=1
εc τc
dce = f (c, ce ), dt
j=1
(12)
where i = 1, . . . , N enumerates the synaptic terminals, and j = 1, . . . , M counts the neurons that are closely located to the glion (generally, N = M). In Eq. (12), c can represent not only single-cell calcium dynamics but also its variation in small networks of glions coupled via gap junctions. Acknowledgments This work was supported by the European Union through the Network of Excellence BioSim (Contract No. LSHB-CT-2004-005137) and RFBR grant 04-0216769. The authors acknowledge fruitful discussions with N. Brazhe, A. Brazhe, and G. Maksimov (Moscow State University).
D.E. Postnov et al. / BioSystems 89 (2007) 84–91
References Araquo, A., Parpura, V., Zanzgiri, R.P., Haydon, P.H., 1999. Tripartite synapses: glia, the unacknowledged. Trends Neurosci. 22, 208. Bonvento, G., Giaume, C., Lorenceau, J., 2002. Neuron–glia interactions: from physiology to behavior. J. Physiol. 96, 167– 168. Carmignoto, G., 2000. Reciprocal communication systems between astrocytes and neurons. Prog. Neurobiol. 62, 561–581. Deitmer, J.W., Rose, C.R., Munch, T., Schmidt, J., Nett, W., Schneider, N.-P., Lohr, C., 1999. Leech glial cell: functional role in a simple nervous system. Glia 28, 175–182. Fields, R.D., Stevens-Graham, B., 2002. New insights into neuron–glia communication. Science 298, 556–562.
91
FitzHugh, R.A., 1961. Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–446. Izhikevich, E.M., 2001. Resonate-and-fire neurons. Neural Netw. J. 14, 883–894. Keener, J., Sneyd, J., 1998. Mathematical Physiology. Springer, New York. Kopell, N., Ermentrout, G.B., Whittington, M.A., Traub, R.D., 2000. Gamma rhythms and beta rhythms have different synchronization properties. Proc. Natl. Acad. Sci. U.S.A. 97, 1867–1872. Nadkarni, S., Jung, P., 2003. Spontaneous oscillations of dressed neurons: a new mechanism for epilepsy?. Phys. Rev. Lett. 91, 268101. Nadkarni, S., Jung, P, 2004. Dressed neurons: modeling neural–glia interactions. Phys. Biol. 1, 35–41. Rojtbak, A.I., 1993. Glia. Nauka, Sankt-Petersburg (in Russian).