- Email: [email protected]
Simulation of CA3 region of hippocampus by kinetic models
BioSystems 48 (1998) 247 – 254
Simulation of CA3 region of hippocampus by kinetic models Franscesco Ventriglia * Istituto di Cibernetica, CNR Via Toiano 6, 80072 Arco Felice (NA), Italy
Abstract The global activities of the hippocampal CA3 field which are induced by stimuli coming from entorhinal cortex have been investigated through computional experiments. The simulations are based on a mathematical model formulated within the frame of a kinetic theory of neural systems. To obtain a description closer to the real biological system, the kinetic model has been supplemented with elements describing the dynamics of neurotransmitters in synaptic clefts. Because high frequency activities occur in the CA3 field, elements for the description of receptor saturation and vesicle depletion phenomena have also been taken into account. The computational experiments show that CA3 model reacts to simulated stimuli from entorhinal sources activating the pyramidal neurons. This activity self-organizes within some limited areas, and spreads to the whole CA3 as a patterned reverberation of short period which requires both excitatory and inhibitory neural populations. The reverberating activity is then stopped by massive inhibition produced by slow-inhibitory neural population. The relationships between the parameters linked to the synaptic dynamic and the elements characterizing these patterns have been investigated. © 1998 Elsevier Science Ireland Ltd. All rights reserved. Keywords: CA3; Kinetic theory; Computer simulation; Reverberation; Memory
1. Introduction The involvement of hippocampus and enthorinal cortex in memory processes is well known though the modalities of their participation are not yet known. In fact, the comprehension of activities and functions of these complex brain structures requires a knowledge ranging from single neurons to entire neural systems. The latter * E-mail: [email protected]
level, which in the past was quite beyond the reach of any experimental apparatus, is now investigated by techniques of real-time imaging, and some important experimental results have been recently reported about the entorhinal– hippocampal system (Iijima et al., 1996). These results, though interesting, fail to provide a clear evidence about the role of single sub-fields of the hippocampus in information processing. Mathematical modeling and computer simulations, based on realistic models of neural structures, can
0303-2647/98/$ - see front matter © 1998 Elsevier Science Ireland Ltd. All rights reserved. PII S0303-2647(98)00071-9
248
F. Ventriglia / BioSystems 48 (1998) 247–254
help both in obtaining some of the information lacking and to stimulate new experimental investigations. A model of the whole CA3 containing elements for a simplified description of the synaptic activity has been formulated by the present author and results of preliminary computational ‘experiments’ were reported in a recent article (Ventriglia, 1998). The approach used is based on a kinetic theory of neural systems and is sharply different from those of others (Traub et al., 1991; Jefferys and Traub, 1994), who also use realistic model based on biological data. In fact, the kinetic theory provides coarse-grained descriptions of the activity of large neuronal systems and investigates space-time features of the activity of whole systems taking into account the most significant neuronal properties (Ventriglia, 1974, 1994). The methods used by the above authors, vice versa, provide thorough descriptions of single neurons and can be applied to the study of small neural populations, but hardly are applicable to whole systems. Moreover, unlike the model of Wilson and Cowan (1972), the kinetic theory takes into account explicitly the velocity distribution of impulses (spikes) in the neural space. In this way it is possible to consider different propagation velocities, to accomplish the de-coupling between short- and long-range effects, to use specific connection topographies of biological neural systems, and to investigate synaptic dynamics. The results of the previous work of the author on CA3 activities showed that the response of CA3 to simulated dentate gyrus (DG) stimuli was rapidly driven by the activation of the pyramidal neurons (in 10 ms), self-organized within some limited areas and spread to the whole CA3 as a patterned reverberation of short period. The shape of the spatio-temporal patterns of neuronal activity was strictly determined by the position and the time of the stimuli hitting CA3 in the first few ( 10) milliseconds of activity. When engaged in patterned reverberatory activity, CA3 was unable to react with external stimuli. This also occurred during the succeeding period of high inhibition. The duration of this period can reach also half a second, depending on the amount of (slow) inhibition produced within the system. To reduce the possibility of computational artifacts, a
more careful control of the mathematical description has been accomplished whose results are shown in the present paper. In the present work, the modeling of neurotransmitter dynamics is presented in more detail. In addition more elements have also been considered in the construction of the kinetic equations. Neurotransmitter depletion and receptor saturation are also taken into account. Some new computational experiments are carried out to control the reliability of previous results with respect to the modeling improvements and these were confirmed.
2. The kinetic model of CA3 Some essential elements for the comprehension of the kinetic model of CA3 are presented below. See Ventriglia (1994, 1998) for a detailed description of the kinetic theory of neural systems. This theory can be partitioned in two stages: a global one, describing electric activity propagation within the system and a local one, describing synaptic activities. They are cross-linked, because the propagation of the electric activity induces changes in neurotransmitter concentrations and these changes produce modifications in the membrane potential of post-synaptic neurons, which can induce firing and propagation of new electric activity. The global stage takes into account the numerous and different types of neuronal populations of a neural system, their densities, their excitatory/inhibitory characteristics, the short and long-range structure of connections, the various types of external input to the system. The local stage considers diverse types of neurotransmitters released in the synaptic clefts, their different kinetics, their different effects on post-synaptic membrane potentials. The present model of CA3 activity took into account three neural populations: pyramidal cells and slow and fast inhibitory cells (Thalmann and Ayala, 1982; Segal, 1990). Several (six) types of impulses were also considered in the model. Two types were conveyed by fibers coming from external sources: granular layer of DG, pyramidal layers of entorhinal cortex (EC). The other three types, short-range impulses, were produced in CA3. One is related to
F. Ventriglia / BioSystems 48 (1998) 247–254
recurrent collaterals of pyramidal cell axons which remain in the neighborhood of the emitting neuron and the other two are related to inhibitory interneurons. Also impulses with long-range effects were considered. They are related to the recurrent collaterals of pyramidal neurons reaching distant zones. To facilitate the mathematical description, the CA3 field was collapsed into a rectangle containing the neural populations. The activity of this neural system was described by the methods of the neural kinetic theory (Ventriglia, 1974, 1994), which considers a neural system as a physical system composed of two families of interacting particles–neurons and impulses, embedded in a bounded volume of space (reduced here to a bounded surface A of a plane). A vector r is related to the position of impulses and of the neurons and a vector 7 is associated to the impulse velocity. The different classes of impulses and neurons are denoted by indices s and s%, respectively. The neurons are characterized by an unique inner variable e —the sub-threshold excitation, linked to the membrane potential. This variable stores the decaying effects of the inputs received from all the connected neurons. The decay occurs towards a characteristic value of the excitation, called resting level er. The variable e is normalized in the interval (0, 1), the values 0 and 1 correspond to the maximum hyperpolarization level (mhl) and to the firing threshold level, respectively. The different types of neurons emit impulses of different types, when their firing threshold is exceeded and they remain thereafter in a state of refractoriness unable to fire for a characteristic time t depending on the type of neuron. The impulses are considered as particles moving along linear paths. Their total number changes in time due to interaction with neurons which absorb impulses (according to an absorption coefficient s) and emit impulses. During the impulse absorption, a quantal amount of specific neurotransmitter is released in the synaptic cleft producing an increase in the related neurotransmitter concentration. This, in turn, induces a change in the sub-threshold excitation e of absorbing (postsynaptic) neurons consisting of an increase (depo-
249
larization) or a decrease (hyperpolarization) of e, according to the excitatory or inhibitory quality of the absorbed impulses. To simulate the external sources of impulses some appropriate functions Zs (r, 7, t) can be chosen giving the number of impulses of type s produced in (r, dr) for unit of time. The short-range, random impulses are emitted by pyramidal and inhibitory neurons as per an emission coefficient ls(r) and diffuse over the surrounding space based on a pdf f *(7). The excitatory, long-range impulses are emitted in each point of the neural field by pyramidal neurons according to an emission coefficient j expressed as an origin-destination matrix: for each couple of points (r, r%) of the field CA3 j(r, r%)
!
"0 =0
if (r, dr) is coupled to (r%, dr%) otherwise. (1)
The number of distant zones reached by these impulses is related to the recurrent collaterals of neurons in the zone of origin. A randomization algorithm, constrained according to the results of Ishizuka et al. (1990), was used for simulations of the CA3 activity (Ventriglia, 1998). In the region between the origin and the destination, the impulses propagate without absorption and the modulus of the velocity can be distributed according a pdf f %( 7 ). Upon reaching the destination, they diffuse with a velocity distributed as per another pdf f 0(7). The absorption coefficient ss,s% depends both on s and s%. The description of the macroscopic neural activity (for example the local frequency of spikes, the local number of firing neurons and the local field potential) requires the knowledge of the following quantities in the volume element (r, dr) at time t for all the different classes s (type of impulse) and s% (type of neuron): the probable number of impulses fs (r, 7, t) dr d7 whose velocities lie in the element d7 around 7 the probable number of neurons gs%(r, e, t) dr de whose excitation lies in (e, e+ de).
F. Ventriglia / BioSystems 48 (1998) 247–254
250
The value e= 0, mhl state, is a special point for the probability distribution function of subthreshold excitation. In fact a finite amount of probability can be associated to it. This probability is described by a separate function g 0s%(r, t) dr, giving the probable number of neurons in hyperpolarization state. The probable number of neurons in refractory state also is described by a separate function g¯s%(r, t) dr. The functions gs% and g¯s% must obey the following principle of conservation of probabilistic mass cs%(r)=g 0s%(r, t)+
&
1
gs%(r, e, t) de +g¯s%(r, t)
(2)
0
where cs%(r) dr is the probable number of neurons of type s% for unit volume of CA3. The absorption rate of impulses of type s which must be differentiated for each type of neuron s%, due to the different role assumed, is given by
&
hs,s%(r, t)= fs (r, 7, t− t0)ss,s%cs%(r) 6 d7,
(3)
where t0 is the synaptic delay. The description of the neurotransmitter dynamic in the synaptic clefts located in an elementary volume near r is carried out by means of two distributions Cs,s%(r, t)
and
C( s,s%(r, t),
which represent the concentrations of neurotransmitters (released by impulses of type s when absorbed by neurons of type s%) which move freely in synaptic clefts (the first distribution) or are linked to receptors on the post-synaptic membrane (the second one). Thus, C( s,s%(r, t) is also related to the concentration of active channels. The equations in Ventriglia (1998) assumed the sudden elimination of the neurotransmitters from the synaptic clefts when they leave the receptorbinding sites. The new equations are based on a two step model: the neurotransmitters leaving the receptors come back in the free moving neurotransmitter pool within the synaptic clefts and they can bind again to the receptors; an active, re-uptake mechanism eliminates neurotransmitters from the synaptic clefts, transporting them (as a whole or as fragments) to pre-synaptic neurons.
These equations are: (Cs,s% (r, t) (t = − kdCs,s%(r, t)−krs,s %Cs,s%(r, t)−kls,s %Cs,s%(r, t) +qas,s %hs,s%(r, t)+kis,s %C( s,s%(r, t)
(4)
(C( s,s% (r, t)= − kis,s %C( s,s%(r, t)+kls,s %Cs,s%(r, t− t( s) (t (5) where t( s is the PSP onset time (we assume that it is different from zero only for slow-IPSPs), kd is the rate constant of neurotransmitters diffusing out of synaptic clefts which will not bind to the post-synaptic membrane (and so they will have no effect on the change of potential), kls,s % is the rate constant of neurotransmitters binding to the postsynaptic membrane receptors which become responsible for the induction of a variation in the inner neuronic excitation through the modification of the opening time of related channels, krs,s % is the rate constant of neurotransmitter re-uptake from synaptic clefts by presynaptic neurons, qas,s % is related to the quantal contribution to the concentration growth of freely moving neurotransmitters in synaptic clefts due to a single impulse absorption, kis,s % is the deactivation constant, i.e. the measure of the rate with which the neurotransmitters lose the control of the mechanisms responsible for the membrane potential change. Because synaptic vesicles (hence, neurotransmitters) and receptors located in r are finite in number, a depletion mechanism, with metabolic recovering, which forbids the unlimited increasing of neurotransmitter concentration and a saturation mechanism for activated receptors were introduced. This is necessary whenever high frequencies of firing are produced within the system (as it occurs in presence of sharp waves). For this purpose we have introduced functions R CV s,s%(r, t) and C s,s%(r, t), which are related to, respectively, the neurotrasmitter concentration in the pre-synaptic vesicle pool (which can be depleted by continuous firing) and the receptor concentration in the post-synaptic neuron pool located in r. The two mechanisms are described by the following equations:
F. Ventriglia / BioSystems 48 (1998) 247–254
(C V s,s% (r, t) (t
251
ets%(r)= % ksC( s,s%(r, t)− % ksC( s,s%(r, t) sI e
= −qas,s %hs,s%(r, t) + kres,s %(1−1/(1 + (exp(−((C V s,s%(r, t)− CVs,s%) − r2)/r1))))
(6)
(C R s,s% (r, t)= − kls,s %Cs,s%(r, t) + kis,s % C( s,s%(r, t) (t
(7)
R where: C V s,s%(r, 0) = CVs,s% and C s,s%(r, 0) = CRs,s%, kres,s% is the rate constant of neurotransmitter forming through metabolic mechanisms, and r1 and r2 are constant. When the function C V s,s%(r, t) reaches the null value (the pre-synaptic vesicles are depleted), no new neurotransmitter can be poured into the synaptic clefts in r and the pre-synaptic firing is without effects. When the function C R s,s%(r, t) is null (the post-synaptic receptors are saturated), no neurotrasmitter can bind to receptors in r. The law of transformation of the concentration variation of membrane-bound neurotransmitters to the inner excitation variation (post-synaptic excitation) was assumed linear in the interval (0, 1). The variation of type s% neuron excitation due to the impulses of type s in an elemental volume (r, dr) at time t is given by
which represents the contribution to the membrane potential variation of membrane-bound neurotransmitter concentration at that time: ens%(r)= ers% + ets%(r). In the case of maximum hyperpolarization, the sub-threshold excitation also reaches the value ens%(r). The probable number of firing neurons in (r, dr) was derived in Ventriglia (1998)
&
Ns%(r, t) 1
=
gs%(r, e, t) de
R(1 − R(DEs %))
+u(ens% − 1)Ns%(r, t− ts%).
(11)
The functions R and u are:
! !
R(x)=
x 0
if x\0 otherwise
(12)
u(x) =
1 0
if x] 0 otherwise
(13)
The following equation, governing the probable number of neurons in mhl state, was also derived: g 0s%(r, t+ D)
&
R( − DEs %)
=
gs%(r, e, t) de
0
DEs,s%(r, t)= ksDC( s,s%(r, t)
(8)
where ks is a constant. Appropriate choices of the parameters in the above equations, producing different PSPs closely related to the PSPs recorded in neurons of CA3 (Miles, (1990)) were made for simulations. The sum of all the contributions both excitatory and inhibitory to the variation of type s% neuron excitation is therefore given by DEs%(r, t)= % DEs,s%(r, t) − % DEs,s%(r, t) sI e
(10)
sI i
(9)
sI i
where Ie and Ii are sets of indices of excitatory and inhibitory impulses in which the whole class of impulses is split. When neurons of type s% exit from refractory state they again are able to react to synaptic effects and the internal excitation reaches a value ens%(r) equal to the sum of resting value e %rs and
+u(− DEs%){g 0s%(r, t)+u(− ens%)Ns%(r, t− ts%)}. (14) The following set of finite-difference kinetic equations was also constructed for the functions fs (r, 7, t), gs%(r, e, t), g¯s%(r, t): fs (r+7D, 7, t+ D) = fs (r, 7, t)− 7 fs (r, 7, t)% cs%(r)ss,s%D s%
+ ls (r)f*(7)Ns (r, t)
&
&
js (r, r%) dr% f %(7%)Ns (r%, t−
+ f 0(7)
A
r−r% ) 7%
d7%+Zs (r, 7, t)D
(15)
gs%(r, e, t+ D) = d(e− ens%)Ns%(r, t− ts%)+ gs%(r, e −DEs%, t) + u(DEs%)d(e− ens%)g 0s%(r, t)
(16)
F. Ventriglia / BioSystems 48 (1998) 247–254
252
g¯s%(r, t +D) = g¯s%(r, t) −Ns%(r, t −ts%) +Ns%(r, t) (17) where d(x) =
!
1 0
if x=0 otherwise.
(18)
This mathematical model has been utilized to carry out a series of computational experiments. The space-time course of some macroscopic parameters (local frequency of spikes, local mean sub-threshold excitation, number of firing neurons), which have close analogy with the in vivo recorded activity of the CA3 field (population spike trains, local field potentials) was analyzed to obtain information on its ability to simulate real hippocampal activity. The meaning of all the parameters was fully explained in Ventriglia (1998) where also the details of the computer simulation setup can be found.
3. Computational experiments and results Some series of computational experiments were carried out where the reaction of the CA3 model to simulated input from mossy fibers of granular neurons of DG and from pyramidal neurons of EC was investigated. In these simulations input from DG was conveyed by layered mossy fibers which stimulate, at different times, different pyramidal strips of CA3, each 3 mm long and 50 mm large. The input to each strip conveyed random volleys whose arrival times were distributed according to a Poisson function. The duration and the amplitude of each volley were Gaussian distributed. Numerical methods discussed in Ventriglia (1988), based on techniques of dynamic programming (Angel and Bellman, 1972) were utilized to solve the kinetic equations quoted in the previous section. Moreover, almost all the parameters maintained the values quoted in Ventriglia (1998), except the description of synaptic dynamics. The changes concerned the presence of new equations for depletion and saturation effects, of a re-uptake term in the synaptic equations, and of variations in the values of the binding coefficient kl as well as of the inactivating
coefficient ki. The other important change, a technical one, was the implementation of a parallelized version of the computer program in a parallel computer, based on two bi-processor nodes. The greater computational power allowed a finer space–time description of CA3 activities. The space and time steps were 50 mm and 0.125 ms, respectively (vs 100 mm and 0.25 ms in Ventriglia, 1998). Three different types of computational experiments were carried out which simulated 750 ms of CA3 activity. For the lack of space, only results for some of the macroscopic variables and for a very small number of times are shown in Fig. 1, where the space behavior of the density of three different types of impulses are represented by three horizontal strips. The impulses originate in CA3 from (going from the bottom): (1) fast inhibitory neurons, (2) slow inhibitory neurons, (3) (short-range axons of) pyramidal neurons. (i) A reduced model, equivalent to the basic model of Ventriglia (1998), was subjected to a first series of simulations for comparisons. In this case the synaptic equations lacked the re-uptake term (krs,s % = 0) and the neurotransmitters leaving the receptor binding sites were instantaneously and completely eliminated from synaptic clefts (i.e. the term kis,s %C( s,s%(r, t) was lacking in Eq. (4)). The saturation and depletion phenomena were not described. In part A of Fig. 1, the global space– time behavior of the CA3 reaction to DG and EC input is shown at four different times. The behavior is similar to the analogous case of Ventriglia (1998) quoted in the Section 1, except the finer resolution. (ii) In a second series of computational experiments the complete set of synaptic Eqs. (4)–(7) was taken into account. In these experiments the value given to the re-uptake coefficient krs,s % was half of that of the binding coefficient kls,s % for all the impulse–neuron couples. In this way, for each couple of neurotransmitter molecules which bound to post-synaptic receptors, one neurotransmitter molecule was transported back to pre-synaptic neurons. Also the saturation effects were active. The part B of Fig. 1 shows the global reaction of the simulated CA3 where the behavior is similar to previous case. About the same time is necessary to induce the reaction, but greater duration and variability are shown. (iii) A greater
F. Ventriglia / BioSystems 48 (1998) 247–254
re-uptake coefficient also was used in a third experiment where krs,s % =kls,s %. Hence, the amount of neurotransmitters bound to post-synaptic receptors is the same as that which is cleared from the synaptic cleft. The results, shown in part C of Fig. 1, demonstrate the induction of a global reaction in CA3. The reaction was delayed, shorter and not well shaped, but involved the whole CA3 and was followed by the same duration of inhibition of the entire field as in case (i) and (ii). The results of the present and of the previous work (Ventriglia, 1998) seem to outline a possible role of CA3 in memory processes. In fact, each series of input from entorhinal cortex and from polymodal cortical areas, having enough power to induce global ‘resonances’ in CA3, are able not only to influence the next stages of the memory process but also to forbid the subsequent information arriving at CA3 (from the same cortical regions) in producing memory effects for a period
253
lasting up to about half a second. This quasi automatic action of resonance and filtering by CA3 requires that the activity of its inhibitory populations and its excitatory input from DG are carefully controlled. In the hippocampus the control on the DG input is exerted by Hilus which is, in turn, controlled by several brain nuclei such as locus coeruleus, diagonal band, raphe, etc. conveying information of high emotional contents. The same nuclei control also the inhibitory populations of CA3. All these nuclei are linked to pain, stress, pleasantness, novelty and other contents which are very important for individual and species survival. In conclusion, the computational results of the simulated model of CA3 and the structural information about the entorhinal– hippocampal system seem to suggest that only cortical activities correlated with sub-cortical activities with important emotional contents have the possibility to save traces in the memory system.
Fig. 1. Space – time behavior of CA3 reactions to simulated input from DG and EC. The three sub-sections (A, B, C) show the results of three different computational experiments discussed in the text. The densities of impulses moving around on the surface of the simulated CA3 are shown in the elementary rectangles. Eighty levels of grey are used, white being the maximum value.
254
F. Ventriglia / BioSystems 48 (1998) 247–254 Thalmann, R.H., Ayala, G.F., 1982. A late increase in potassium conductance follows synaptic stimulation of granule neurons of the dentate gyrus. Neurosci. Lett. 29, 243 – 248. Traub, R.D., Wong, R.K., Miles, R.S., Michelson, H., 1991. A model of a CA3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances. J. Neurophysiol. 66, 635 – 650. Ventriglia, F., 1974. Kinetic approach to neural systems. I. Bull. Math. Biol. 36, 534 – 544. Ventriglia, F., 1988. Computational simulation of activity of cortical-like neural systems. Bull. Math. Biol. 50, 143 – 185. Ventriglia, F., 1994. Towards a kinetic theory of cortical-like neural fields. In: Ventriglia, F. (Ed.), Neural Modeling and Neural Networks. Pergamon Press, Oxford, pp. 217 – 249. Ventriglia, F., 1998. Computational experiments support a competitive function in the CA3 region of the hippocampus. Bull. Math. Biol. 60, 373 – 407. Wilson, H.R., Cowan, J.D., 1972. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1 – 24.
References Angel, E., Bellman, R., 1972. Dynamic programming and partial differential equations. Academic Press, New York. Iijima, T., Witter, M.P., Ichikawa, M., Tominaga, T., Kajiwara, R., Matsumoto, G., 1996. Entorhinal–hippocampal interactions revealed by real-time imaging. Science 272, 1176–1179. Ishizuka, N., Weber, J., Amaral, D.G., 1990. Organization of intra-hippocampal projections originating from CA3 pyramidal cells in the rat. J. Comp. Neurol. 295, 580–623. Jefferys, J.G.R., Traub, R.D., 1994. Mechanisms responsible for epilepsy in hippocampal slices predispose the brain to collective oscillations. In: Ventriglia, F. (Ed.), Neural Modeling and Neural Networks. Pergamon Press, Oxford, pp. 111 – 127. Miles, R., 1990. Synaptic excitation of inhibitory cells by single CA3 hippocampal pyramidal cells of guinea-pig hippocampus in vitro. J. Physiol. 431, 659–676. Segal, M., 1990. A subset of local interneurons generate slow inhibitory postsynaptic potentials in hippocampal neurons. Brain Res. 511, 163 –164.
.
Simulation of CA3 region of hippocampus by kinetic models Franscesco Ventriglia * Istituto di Cibernetica, CNR Via Toiano 6, 80072 Arco Felice (NA), Italy
Abstract The global activities of the hippocampal CA3 field which are induced by stimuli coming from entorhinal cortex have been investigated through computional experiments. The simulations are based on a mathematical model formulated within the frame of a kinetic theory of neural systems. To obtain a description closer to the real biological system, the kinetic model has been supplemented with elements describing the dynamics of neurotransmitters in synaptic clefts. Because high frequency activities occur in the CA3 field, elements for the description of receptor saturation and vesicle depletion phenomena have also been taken into account. The computational experiments show that CA3 model reacts to simulated stimuli from entorhinal sources activating the pyramidal neurons. This activity self-organizes within some limited areas, and spreads to the whole CA3 as a patterned reverberation of short period which requires both excitatory and inhibitory neural populations. The reverberating activity is then stopped by massive inhibition produced by slow-inhibitory neural population. The relationships between the parameters linked to the synaptic dynamic and the elements characterizing these patterns have been investigated. © 1998 Elsevier Science Ireland Ltd. All rights reserved. Keywords: CA3; Kinetic theory; Computer simulation; Reverberation; Memory
1. Introduction The involvement of hippocampus and enthorinal cortex in memory processes is well known though the modalities of their participation are not yet known. In fact, the comprehension of activities and functions of these complex brain structures requires a knowledge ranging from single neurons to entire neural systems. The latter * E-mail: [email protected]
level, which in the past was quite beyond the reach of any experimental apparatus, is now investigated by techniques of real-time imaging, and some important experimental results have been recently reported about the entorhinal– hippocampal system (Iijima et al., 1996). These results, though interesting, fail to provide a clear evidence about the role of single sub-fields of the hippocampus in information processing. Mathematical modeling and computer simulations, based on realistic models of neural structures, can
0303-2647/98/$ - see front matter © 1998 Elsevier Science Ireland Ltd. All rights reserved. PII S0303-2647(98)00071-9
248
F. Ventriglia / BioSystems 48 (1998) 247–254
help both in obtaining some of the information lacking and to stimulate new experimental investigations. A model of the whole CA3 containing elements for a simplified description of the synaptic activity has been formulated by the present author and results of preliminary computational ‘experiments’ were reported in a recent article (Ventriglia, 1998). The approach used is based on a kinetic theory of neural systems and is sharply different from those of others (Traub et al., 1991; Jefferys and Traub, 1994), who also use realistic model based on biological data. In fact, the kinetic theory provides coarse-grained descriptions of the activity of large neuronal systems and investigates space-time features of the activity of whole systems taking into account the most significant neuronal properties (Ventriglia, 1974, 1994). The methods used by the above authors, vice versa, provide thorough descriptions of single neurons and can be applied to the study of small neural populations, but hardly are applicable to whole systems. Moreover, unlike the model of Wilson and Cowan (1972), the kinetic theory takes into account explicitly the velocity distribution of impulses (spikes) in the neural space. In this way it is possible to consider different propagation velocities, to accomplish the de-coupling between short- and long-range effects, to use specific connection topographies of biological neural systems, and to investigate synaptic dynamics. The results of the previous work of the author on CA3 activities showed that the response of CA3 to simulated dentate gyrus (DG) stimuli was rapidly driven by the activation of the pyramidal neurons (in 10 ms), self-organized within some limited areas and spread to the whole CA3 as a patterned reverberation of short period. The shape of the spatio-temporal patterns of neuronal activity was strictly determined by the position and the time of the stimuli hitting CA3 in the first few ( 10) milliseconds of activity. When engaged in patterned reverberatory activity, CA3 was unable to react with external stimuli. This also occurred during the succeeding period of high inhibition. The duration of this period can reach also half a second, depending on the amount of (slow) inhibition produced within the system. To reduce the possibility of computational artifacts, a
more careful control of the mathematical description has been accomplished whose results are shown in the present paper. In the present work, the modeling of neurotransmitter dynamics is presented in more detail. In addition more elements have also been considered in the construction of the kinetic equations. Neurotransmitter depletion and receptor saturation are also taken into account. Some new computational experiments are carried out to control the reliability of previous results with respect to the modeling improvements and these were confirmed.
2. The kinetic model of CA3 Some essential elements for the comprehension of the kinetic model of CA3 are presented below. See Ventriglia (1994, 1998) for a detailed description of the kinetic theory of neural systems. This theory can be partitioned in two stages: a global one, describing electric activity propagation within the system and a local one, describing synaptic activities. They are cross-linked, because the propagation of the electric activity induces changes in neurotransmitter concentrations and these changes produce modifications in the membrane potential of post-synaptic neurons, which can induce firing and propagation of new electric activity. The global stage takes into account the numerous and different types of neuronal populations of a neural system, their densities, their excitatory/inhibitory characteristics, the short and long-range structure of connections, the various types of external input to the system. The local stage considers diverse types of neurotransmitters released in the synaptic clefts, their different kinetics, their different effects on post-synaptic membrane potentials. The present model of CA3 activity took into account three neural populations: pyramidal cells and slow and fast inhibitory cells (Thalmann and Ayala, 1982; Segal, 1990). Several (six) types of impulses were also considered in the model. Two types were conveyed by fibers coming from external sources: granular layer of DG, pyramidal layers of entorhinal cortex (EC). The other three types, short-range impulses, were produced in CA3. One is related to
F. Ventriglia / BioSystems 48 (1998) 247–254
recurrent collaterals of pyramidal cell axons which remain in the neighborhood of the emitting neuron and the other two are related to inhibitory interneurons. Also impulses with long-range effects were considered. They are related to the recurrent collaterals of pyramidal neurons reaching distant zones. To facilitate the mathematical description, the CA3 field was collapsed into a rectangle containing the neural populations. The activity of this neural system was described by the methods of the neural kinetic theory (Ventriglia, 1974, 1994), which considers a neural system as a physical system composed of two families of interacting particles–neurons and impulses, embedded in a bounded volume of space (reduced here to a bounded surface A of a plane). A vector r is related to the position of impulses and of the neurons and a vector 7 is associated to the impulse velocity. The different classes of impulses and neurons are denoted by indices s and s%, respectively. The neurons are characterized by an unique inner variable e —the sub-threshold excitation, linked to the membrane potential. This variable stores the decaying effects of the inputs received from all the connected neurons. The decay occurs towards a characteristic value of the excitation, called resting level er. The variable e is normalized in the interval (0, 1), the values 0 and 1 correspond to the maximum hyperpolarization level (mhl) and to the firing threshold level, respectively. The different types of neurons emit impulses of different types, when their firing threshold is exceeded and they remain thereafter in a state of refractoriness unable to fire for a characteristic time t depending on the type of neuron. The impulses are considered as particles moving along linear paths. Their total number changes in time due to interaction with neurons which absorb impulses (according to an absorption coefficient s) and emit impulses. During the impulse absorption, a quantal amount of specific neurotransmitter is released in the synaptic cleft producing an increase in the related neurotransmitter concentration. This, in turn, induces a change in the sub-threshold excitation e of absorbing (postsynaptic) neurons consisting of an increase (depo-
249
larization) or a decrease (hyperpolarization) of e, according to the excitatory or inhibitory quality of the absorbed impulses. To simulate the external sources of impulses some appropriate functions Zs (r, 7, t) can be chosen giving the number of impulses of type s produced in (r, dr) for unit of time. The short-range, random impulses are emitted by pyramidal and inhibitory neurons as per an emission coefficient ls(r) and diffuse over the surrounding space based on a pdf f *(7). The excitatory, long-range impulses are emitted in each point of the neural field by pyramidal neurons according to an emission coefficient j expressed as an origin-destination matrix: for each couple of points (r, r%) of the field CA3 j(r, r%)
!
"0 =0
if (r, dr) is coupled to (r%, dr%) otherwise. (1)
The number of distant zones reached by these impulses is related to the recurrent collaterals of neurons in the zone of origin. A randomization algorithm, constrained according to the results of Ishizuka et al. (1990), was used for simulations of the CA3 activity (Ventriglia, 1998). In the region between the origin and the destination, the impulses propagate without absorption and the modulus of the velocity can be distributed according a pdf f %( 7 ). Upon reaching the destination, they diffuse with a velocity distributed as per another pdf f 0(7). The absorption coefficient ss,s% depends both on s and s%. The description of the macroscopic neural activity (for example the local frequency of spikes, the local number of firing neurons and the local field potential) requires the knowledge of the following quantities in the volume element (r, dr) at time t for all the different classes s (type of impulse) and s% (type of neuron): the probable number of impulses fs (r, 7, t) dr d7 whose velocities lie in the element d7 around 7 the probable number of neurons gs%(r, e, t) dr de whose excitation lies in (e, e+ de).
F. Ventriglia / BioSystems 48 (1998) 247–254
250
The value e= 0, mhl state, is a special point for the probability distribution function of subthreshold excitation. In fact a finite amount of probability can be associated to it. This probability is described by a separate function g 0s%(r, t) dr, giving the probable number of neurons in hyperpolarization state. The probable number of neurons in refractory state also is described by a separate function g¯s%(r, t) dr. The functions gs% and g¯s% must obey the following principle of conservation of probabilistic mass cs%(r)=g 0s%(r, t)+
&
1
gs%(r, e, t) de +g¯s%(r, t)
(2)
0
where cs%(r) dr is the probable number of neurons of type s% for unit volume of CA3. The absorption rate of impulses of type s which must be differentiated for each type of neuron s%, due to the different role assumed, is given by
&
hs,s%(r, t)= fs (r, 7, t− t0)ss,s%cs%(r) 6 d7,
(3)
where t0 is the synaptic delay. The description of the neurotransmitter dynamic in the synaptic clefts located in an elementary volume near r is carried out by means of two distributions Cs,s%(r, t)
and
C( s,s%(r, t),
which represent the concentrations of neurotransmitters (released by impulses of type s when absorbed by neurons of type s%) which move freely in synaptic clefts (the first distribution) or are linked to receptors on the post-synaptic membrane (the second one). Thus, C( s,s%(r, t) is also related to the concentration of active channels. The equations in Ventriglia (1998) assumed the sudden elimination of the neurotransmitters from the synaptic clefts when they leave the receptorbinding sites. The new equations are based on a two step model: the neurotransmitters leaving the receptors come back in the free moving neurotransmitter pool within the synaptic clefts and they can bind again to the receptors; an active, re-uptake mechanism eliminates neurotransmitters from the synaptic clefts, transporting them (as a whole or as fragments) to pre-synaptic neurons.
These equations are: (Cs,s% (r, t) (t = − kdCs,s%(r, t)−krs,s %Cs,s%(r, t)−kls,s %Cs,s%(r, t) +qas,s %hs,s%(r, t)+kis,s %C( s,s%(r, t)
(4)
(C( s,s% (r, t)= − kis,s %C( s,s%(r, t)+kls,s %Cs,s%(r, t− t( s) (t (5) where t( s is the PSP onset time (we assume that it is different from zero only for slow-IPSPs), kd is the rate constant of neurotransmitters diffusing out of synaptic clefts which will not bind to the post-synaptic membrane (and so they will have no effect on the change of potential), kls,s % is the rate constant of neurotransmitters binding to the postsynaptic membrane receptors which become responsible for the induction of a variation in the inner neuronic excitation through the modification of the opening time of related channels, krs,s % is the rate constant of neurotransmitter re-uptake from synaptic clefts by presynaptic neurons, qas,s % is related to the quantal contribution to the concentration growth of freely moving neurotransmitters in synaptic clefts due to a single impulse absorption, kis,s % is the deactivation constant, i.e. the measure of the rate with which the neurotransmitters lose the control of the mechanisms responsible for the membrane potential change. Because synaptic vesicles (hence, neurotransmitters) and receptors located in r are finite in number, a depletion mechanism, with metabolic recovering, which forbids the unlimited increasing of neurotransmitter concentration and a saturation mechanism for activated receptors were introduced. This is necessary whenever high frequencies of firing are produced within the system (as it occurs in presence of sharp waves). For this purpose we have introduced functions R CV s,s%(r, t) and C s,s%(r, t), which are related to, respectively, the neurotrasmitter concentration in the pre-synaptic vesicle pool (which can be depleted by continuous firing) and the receptor concentration in the post-synaptic neuron pool located in r. The two mechanisms are described by the following equations:
F. Ventriglia / BioSystems 48 (1998) 247–254
(C V s,s% (r, t) (t
251
ets%(r)= % ksC( s,s%(r, t)− % ksC( s,s%(r, t) sI e
= −qas,s %hs,s%(r, t) + kres,s %(1−1/(1 + (exp(−((C V s,s%(r, t)− CVs,s%) − r2)/r1))))
(6)
(C R s,s% (r, t)= − kls,s %Cs,s%(r, t) + kis,s % C( s,s%(r, t) (t
(7)
R where: C V s,s%(r, 0) = CVs,s% and C s,s%(r, 0) = CRs,s%, kres,s% is the rate constant of neurotransmitter forming through metabolic mechanisms, and r1 and r2 are constant. When the function C V s,s%(r, t) reaches the null value (the pre-synaptic vesicles are depleted), no new neurotransmitter can be poured into the synaptic clefts in r and the pre-synaptic firing is without effects. When the function C R s,s%(r, t) is null (the post-synaptic receptors are saturated), no neurotrasmitter can bind to receptors in r. The law of transformation of the concentration variation of membrane-bound neurotransmitters to the inner excitation variation (post-synaptic excitation) was assumed linear in the interval (0, 1). The variation of type s% neuron excitation due to the impulses of type s in an elemental volume (r, dr) at time t is given by
which represents the contribution to the membrane potential variation of membrane-bound neurotransmitter concentration at that time: ens%(r)= ers% + ets%(r). In the case of maximum hyperpolarization, the sub-threshold excitation also reaches the value ens%(r). The probable number of firing neurons in (r, dr) was derived in Ventriglia (1998)
&
Ns%(r, t) 1
=
gs%(r, e, t) de
R(1 − R(DEs %))
+u(ens% − 1)Ns%(r, t− ts%).
(11)
The functions R and u are:
! !
R(x)=
x 0
if x\0 otherwise
(12)
u(x) =
1 0
if x] 0 otherwise
(13)
The following equation, governing the probable number of neurons in mhl state, was also derived: g 0s%(r, t+ D)
&
R( − DEs %)
=
gs%(r, e, t) de
0
DEs,s%(r, t)= ksDC( s,s%(r, t)
(8)
where ks is a constant. Appropriate choices of the parameters in the above equations, producing different PSPs closely related to the PSPs recorded in neurons of CA3 (Miles, (1990)) were made for simulations. The sum of all the contributions both excitatory and inhibitory to the variation of type s% neuron excitation is therefore given by DEs%(r, t)= % DEs,s%(r, t) − % DEs,s%(r, t) sI e
(10)
sI i
(9)
sI i
where Ie and Ii are sets of indices of excitatory and inhibitory impulses in which the whole class of impulses is split. When neurons of type s% exit from refractory state they again are able to react to synaptic effects and the internal excitation reaches a value ens%(r) equal to the sum of resting value e %rs and
+u(− DEs%){g 0s%(r, t)+u(− ens%)Ns%(r, t− ts%)}. (14) The following set of finite-difference kinetic equations was also constructed for the functions fs (r, 7, t), gs%(r, e, t), g¯s%(r, t): fs (r+7D, 7, t+ D) = fs (r, 7, t)− 7 fs (r, 7, t)% cs%(r)ss,s%D s%
+ ls (r)f*(7)Ns (r, t)
&
&
js (r, r%) dr% f %(7%)Ns (r%, t−
+ f 0(7)
A
r−r% ) 7%
d7%+Zs (r, 7, t)D
(15)
gs%(r, e, t+ D) = d(e− ens%)Ns%(r, t− ts%)+ gs%(r, e −DEs%, t) + u(DEs%)d(e− ens%)g 0s%(r, t)
(16)
F. Ventriglia / BioSystems 48 (1998) 247–254
252
g¯s%(r, t +D) = g¯s%(r, t) −Ns%(r, t −ts%) +Ns%(r, t) (17) where d(x) =
!
1 0
if x=0 otherwise.
(18)
This mathematical model has been utilized to carry out a series of computational experiments. The space-time course of some macroscopic parameters (local frequency of spikes, local mean sub-threshold excitation, number of firing neurons), which have close analogy with the in vivo recorded activity of the CA3 field (population spike trains, local field potentials) was analyzed to obtain information on its ability to simulate real hippocampal activity. The meaning of all the parameters was fully explained in Ventriglia (1998) where also the details of the computer simulation setup can be found.
3. Computational experiments and results Some series of computational experiments were carried out where the reaction of the CA3 model to simulated input from mossy fibers of granular neurons of DG and from pyramidal neurons of EC was investigated. In these simulations input from DG was conveyed by layered mossy fibers which stimulate, at different times, different pyramidal strips of CA3, each 3 mm long and 50 mm large. The input to each strip conveyed random volleys whose arrival times were distributed according to a Poisson function. The duration and the amplitude of each volley were Gaussian distributed. Numerical methods discussed in Ventriglia (1988), based on techniques of dynamic programming (Angel and Bellman, 1972) were utilized to solve the kinetic equations quoted in the previous section. Moreover, almost all the parameters maintained the values quoted in Ventriglia (1998), except the description of synaptic dynamics. The changes concerned the presence of new equations for depletion and saturation effects, of a re-uptake term in the synaptic equations, and of variations in the values of the binding coefficient kl as well as of the inactivating
coefficient ki. The other important change, a technical one, was the implementation of a parallelized version of the computer program in a parallel computer, based on two bi-processor nodes. The greater computational power allowed a finer space–time description of CA3 activities. The space and time steps were 50 mm and 0.125 ms, respectively (vs 100 mm and 0.25 ms in Ventriglia, 1998). Three different types of computational experiments were carried out which simulated 750 ms of CA3 activity. For the lack of space, only results for some of the macroscopic variables and for a very small number of times are shown in Fig. 1, where the space behavior of the density of three different types of impulses are represented by three horizontal strips. The impulses originate in CA3 from (going from the bottom): (1) fast inhibitory neurons, (2) slow inhibitory neurons, (3) (short-range axons of) pyramidal neurons. (i) A reduced model, equivalent to the basic model of Ventriglia (1998), was subjected to a first series of simulations for comparisons. In this case the synaptic equations lacked the re-uptake term (krs,s % = 0) and the neurotransmitters leaving the receptor binding sites were instantaneously and completely eliminated from synaptic clefts (i.e. the term kis,s %C( s,s%(r, t) was lacking in Eq. (4)). The saturation and depletion phenomena were not described. In part A of Fig. 1, the global space– time behavior of the CA3 reaction to DG and EC input is shown at four different times. The behavior is similar to the analogous case of Ventriglia (1998) quoted in the Section 1, except the finer resolution. (ii) In a second series of computational experiments the complete set of synaptic Eqs. (4)–(7) was taken into account. In these experiments the value given to the re-uptake coefficient krs,s % was half of that of the binding coefficient kls,s % for all the impulse–neuron couples. In this way, for each couple of neurotransmitter molecules which bound to post-synaptic receptors, one neurotransmitter molecule was transported back to pre-synaptic neurons. Also the saturation effects were active. The part B of Fig. 1 shows the global reaction of the simulated CA3 where the behavior is similar to previous case. About the same time is necessary to induce the reaction, but greater duration and variability are shown. (iii) A greater
F. Ventriglia / BioSystems 48 (1998) 247–254
re-uptake coefficient also was used in a third experiment where krs,s % =kls,s %. Hence, the amount of neurotransmitters bound to post-synaptic receptors is the same as that which is cleared from the synaptic cleft. The results, shown in part C of Fig. 1, demonstrate the induction of a global reaction in CA3. The reaction was delayed, shorter and not well shaped, but involved the whole CA3 and was followed by the same duration of inhibition of the entire field as in case (i) and (ii). The results of the present and of the previous work (Ventriglia, 1998) seem to outline a possible role of CA3 in memory processes. In fact, each series of input from entorhinal cortex and from polymodal cortical areas, having enough power to induce global ‘resonances’ in CA3, are able not only to influence the next stages of the memory process but also to forbid the subsequent information arriving at CA3 (from the same cortical regions) in producing memory effects for a period
253
lasting up to about half a second. This quasi automatic action of resonance and filtering by CA3 requires that the activity of its inhibitory populations and its excitatory input from DG are carefully controlled. In the hippocampus the control on the DG input is exerted by Hilus which is, in turn, controlled by several brain nuclei such as locus coeruleus, diagonal band, raphe, etc. conveying information of high emotional contents. The same nuclei control also the inhibitory populations of CA3. All these nuclei are linked to pain, stress, pleasantness, novelty and other contents which are very important for individual and species survival. In conclusion, the computational results of the simulated model of CA3 and the structural information about the entorhinal– hippocampal system seem to suggest that only cortical activities correlated with sub-cortical activities with important emotional contents have the possibility to save traces in the memory system.
Fig. 1. Space – time behavior of CA3 reactions to simulated input from DG and EC. The three sub-sections (A, B, C) show the results of three different computational experiments discussed in the text. The densities of impulses moving around on the surface of the simulated CA3 are shown in the elementary rectangles. Eighty levels of grey are used, white being the maximum value.
254
F. Ventriglia / BioSystems 48 (1998) 247–254 Thalmann, R.H., Ayala, G.F., 1982. A late increase in potassium conductance follows synaptic stimulation of granule neurons of the dentate gyrus. Neurosci. Lett. 29, 243 – 248. Traub, R.D., Wong, R.K., Miles, R.S., Michelson, H., 1991. A model of a CA3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances. J. Neurophysiol. 66, 635 – 650. Ventriglia, F., 1974. Kinetic approach to neural systems. I. Bull. Math. Biol. 36, 534 – 544. Ventriglia, F., 1988. Computational simulation of activity of cortical-like neural systems. Bull. Math. Biol. 50, 143 – 185. Ventriglia, F., 1994. Towards a kinetic theory of cortical-like neural fields. In: Ventriglia, F. (Ed.), Neural Modeling and Neural Networks. Pergamon Press, Oxford, pp. 217 – 249. Ventriglia, F., 1998. Computational experiments support a competitive function in the CA3 region of the hippocampus. Bull. Math. Biol. 60, 373 – 407. Wilson, H.R., Cowan, J.D., 1972. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1 – 24.
References Angel, E., Bellman, R., 1972. Dynamic programming and partial differential equations. Academic Press, New York. Iijima, T., Witter, M.P., Ichikawa, M., Tominaga, T., Kajiwara, R., Matsumoto, G., 1996. Entorhinal–hippocampal interactions revealed by real-time imaging. Science 272, 1176–1179. Ishizuka, N., Weber, J., Amaral, D.G., 1990. Organization of intra-hippocampal projections originating from CA3 pyramidal cells in the rat. J. Comp. Neurol. 295, 580–623. Jefferys, J.G.R., Traub, R.D., 1994. Mechanisms responsible for epilepsy in hippocampal slices predispose the brain to collective oscillations. In: Ventriglia, F. (Ed.), Neural Modeling and Neural Networks. Pergamon Press, Oxford, pp. 111 – 127. Miles, R., 1990. Synaptic excitation of inhibitory cells by single CA3 hippocampal pyramidal cells of guinea-pig hippocampus in vitro. J. Physiol. 431, 659–676. Segal, M., 1990. A subset of local interneurons generate slow inhibitory postsynaptic potentials in hippocampal neurons. Brain Res. 511, 163 –164.
.