Dynamical questions in volume transmission

Abstract

In volume transmission (or neuromodulation) neurons do not make one-to-one connections to other neurons, but instead simply release neurotransmitter into the extracellular space from numerous varicosities. Many well-known neurotransmitters including serotonin (5HT), dopamine (DA), histamine (HA), Gamma-Aminobutyric Acid (GABA) and acetylcholine (ACh) participate in volume transmission. Typically, the cell bodies are in one volume and the axons project to a distant volume in the brain releasing the neurotransmitter there. We introduce volume transmission and describe mathematically two natural homeostatic mechanisms. In some brain regions several neurotransmitters in the extracellular space affect each other's release. We investigate the dynamics created by this comodulation in two different cases: serotonin and histamine; and the comodulation of 4 neurotransmitters in the striatum and we compare to experimental data. This kind of comodulation poses new dynamical questions as well as the question of how these biochemical networks influence the electrophysiological networks in the brain.

Keywords:

• Volume transmission
• neuromodulation
• dynamics
• brain
• chemistry
• comodulation

AMS Classifications:

• 34A37
• 37F46
• 92C20
• 92C40

1. Introduction

It is natural for us to think of the brain as a large computational device that processes information analogously to a computer. In this view, the basic elements are the neurons that receive inputs from other neurons and, via action potentials, send information to other neurons. There are then two fundamental classes of biological and mathematical questions. How do individual neurons receive and process their inputs and decide when to fire? How do connected sets of neurons perform information processing functions that individual neurons cannot do? This electrophysiological point of view is natural for several reasons. If brains are like computational devices then we can use computational algorithms as metaphors and examples of what must be going on in the brain. Secondly, the electrophysiological point of view fits well with our modern scientific method of trying to understand complex behaviour in the large as the interaction of many fundamental parts (the neurons) whose behaviour we understand very well. The electrophysiological point of view is also perfect for mathematical analysis and computation. One need not deal with the messy details of cell biology, the existence of more than 50 identified neurotransmitters, or the fact that neurons come in a bewildering variety of different morphological and physiological types. All of these things appear, if they appear at all, as parameters in models of neurons, or as parameters in local or global network simulations.

However, many neurons are not engaged in one-to-one electrophysiological transmission with other neurons [17]. Instead, groups of neurons whose cell bodies are in one region of the brain project densely to a distant region in the brain in which they don't make synapses but, instead, release their neurotransmitter from a plethora of small varicosities into the extracellular space; see Figure 1. The neurotransmitter in the extracellular space binds to receptors on neurons and other cells in the distant region affecting their behaviour. This is called ‘volume transmission’ or ‘neuromodulation’ and it has been known since the 1960s [17], but evidence has been accumulating that it plays a very important role in brain function and control [19]. The importance of volume transmission has been emphasized by Booth [13] and Rubin [31].

Figure 1. Volume transmission. The cell bodies of neurons in the yellow brain region send their axons to the distant blue volume where they don't make synapses but instead release neurotransmitter from numerous small varicosities into the extracellular space. There, the neurotransmitter binds to receptors on neurons and glial cells changing their behaviour. This is called ‘volume transmission’ or ‘neuromodulation.’

For example, almost all of the cell bodies of serotonin(5HT) neurons in the brain are contained in two small nuclei, the Dorsal Raphe Nucleus(DRN) and the Medial Raphe Nucleus(MRN) in the brainstem. The axons of these neurons project to many distant brain regions and many of these regions send projections back to the DRN and the MRN [28]. Most or all of these projections are thought to be volume transmission. For instance, 5HT in the striatum, a part of the basal ganglia, binds to receptors on dopamine(DA) synapses increasing the release of DA per action potential. In fact, there are more than 50 different 5HT receptors in the brain [2], suggesting that this great variety allows the concentration of 5HT to modulate neurons in different ways in different regions depending on which 5HT receptors the cells in the region express. The neurons engaged in volume transmission fire action potentials, of course, but the action potentials do not carry information in the usual sense; instead, their purpose is to allow the neurons to project biochemical changes over long distances in the brain.

It has been usual in the literature to refer to the volume transmission of different neurotransmitters as the ‘serotonergic’ system and the ‘dopaminergic’ system and so forth, as though these systems were separate. However, it has become clear [1, 29, 33, 36] that the situation is even more complicated because the different volume transmission systems affect each other and that they ‘comodulate’ or ‘coregulate’ each other. This raises many new biological and mathematical questions. In Sections 3 and 4 we present new work on this comodulation. In Section 2, we review older work on homeostatic mechanisms in volume transmission that form the basis for the work in Sections 3 and 4. In Section 6, we discuss the mathematical models that we use, so that others can investigate the dynamical questions themselves. Finally, in the Discussion we describe the many new mathematical questions that arise in comodulation.

2. Homeostatic mechanisms in volume transmission

In this section we discuss two homeostatic mechanisms in volume transmission that stabilize neurotransmitter concentrations in target regions, passive stabilization and control by autoreceptors. Passive stabilization refers to the following idea. Suppose that half the cells in the yellow region in Figure 1 die. What will happen to the concentration of the neurotransmitter in the blue projection region? At first one might guess that the concentration would decrease by half, but that is not correct.

Released neurotransmitter is taken back up into the black varicosities by transporters and the observed concentration reflects a balance between release and reuptake. When half the cells in the yellow region die then there will be approximately half as much release but also half as many reuptake transporters, so the balance remains the same and the concentration shouldn't change. That's not quite true because there's another effect. In addition to reuptake into varicosities, neurotransmitter released in the blue region can be transported into glial cells or taken up into blood vessels or simply diffuse out of the tissue. As the remaining varicosities get sparser and sparser (as cells in the yellow region die) the probability of these removal processes goes way up because it is hard for a released neurotransmitter molecule to find a transporter and be taken back up into varicosities. The observed concentration depends on both of these effects and so requires a computation. We use the model for serotonin and histamine neurons described in Section 5. We compute the concentration of 5HT in the extracellular space in the blue region as a function of the fraction, f, of serotonin neurons left alive in the yellow region. The result can be seen in Figure 2. The concentration of 5HT in the extracellular space shows remarkable homeostasis as cells die until about 80% of the cells have died and then the concentration plunges to zero.

Figure 2. Passive Stabilization. As the fraction of cells alive decreases, the concentration of 5HT in the projection region remains quite stable until 80% of the cells have died. See Section 5 for a description of the model.

This passive stabilization is a general phenomenon for volume transmission that we have illustrated using our serotonin model. The idea behind it was first explained in the case of DA by Bergstrom and Garris [5] and supported by computations with our DA model [34]. The DA case is very important, because in Parkinson's disease, DA neurons in the substantia nigra pars compacta (SNc) die. These DA cells project to the striatum in the basal ganglia, which is a motor control region of the brain, where it is known that abnormally low concentrations of DA cause the symptoms of Parkinson's disease. Passive stabilization explains why Parkinson's symptoms are generally not observed until 80% of the cells in the SNc have died. This is supported by a large number of experimental and clinical studies [3, 6, 10, 14, 16, 25, 30, 35, 38, 41]. There is a substantial amount of individual variability and the 80% value seems to fit well with most summary data in the literature.

To explain the role of autoreceptors we need to discuss the synthesis, storage, release and reuptake of monoamine neurotransmitters. We'll use serotonin and our example; see Figure 10 in Section 6. Serotonin (5HT) is synthesized in varicosities from the amino acid tryptophan in the blood, which crosses the blood-brain barrier and enters the extracellular fluid from which it is transported into the varicosity. There, the enzyme tryptophan hydroxylase (TPH) adds an OH group and the enzyme Aromatic Amino acid Decarboxylase (AADC) cuts off the carboxyl group to make cytosolic 5HT. The cytosolic 5HT is packaged into vesicles (v5HT) by the monoamine transporter (VMAT), and when the action potential passes by the varicosity some vesicles move to the boundary and release 5HT into the extracellular space (e5HT). There it diffuses and is also transported back into varicosities by the serotonin transporter (SERT). On the boundary of the varicosity are large protein receptors (5HT1b), called autoreceptors, to which the extracellular 5HT can bind indicating to the cell the concentration in the extracellular space. If the concentration is too high, the autoreceptors inhibit release as indicated in Figure 10. And, when the extracellular 5HT is lower than the steady state concentration, the autoreceptors do not inhibit release as much thus causing 5HT concentrations to rise. This kind of end-product inhibition is very common in biochemical systems. Here it is clearly a mechanism to stabilize the concentration of 5HT in the extracellular space. Figure 3 shows the effect of the serotonin autoreceptor. In the presence of the autoreceptor the perturbation in the extracellular space caused by the release from the 5HT neuron is much smaller and comes back to baseline more quickly.

Figure 3. Autoreceptors. The curves show the concentration of 5HT in the extracellular space in response to a 1.5 s stimulation (black bar) without (blue) and with (red) autoreceptors. See Section 6 for a description of model.

The story we have just described is well-known [17]. However, the detailed mechanisms by which the autoreceptors influence release are still under study. Here we have used a very simple model for the autoreceptor where the concentration of 5HT in the extracellular space directly influences release (see Section 6). In fact, the dynamics of the autoreceptors are more complicated and more interesting. They are G-protein coupled receptors that stimulate cellular machinery. In some of our models [7, 9] we give more detailed models of the G-protein dynamics.

3. Comodulation of serotonin and histamine

In previous work, we developed separate mathematical models for histamine and serotonin synthesis, storage, release, reuptake including an H3 autoreceptor on the histamine varicosity and a 5HT1b serotonin autoreceptor on the serotonin varicosity; see Figure 10. In recent years, it has been shown that there are H3 receptors also on 5HT varicosities and when HA binds it inhibits 5HT release [7, 9, 37]. There is also evidence of a 5HT receptor on HA varicosities that stimulates HA release when 5HT levels are above equilibrium and inhibits HA release when 5HT levels are below equilibrium in the extracellular space [15, 27, 29]. This means that these two neurotransmitters influence each other's release in areas of the brain where both types of varicosities are present. Not only does extracellular HA inhibit 5HT release, but extracellular 5HT stimulates HA release. This leads us to the following question: If either neuron is stimulated, what are the dynamics of the two neurotransmitters in the extracellular space?

To answer this question, we combined our previous 5HT and HA programmes that include the HA receptor on the 5HT varicosity and added the 5HT receptor on the HA varicosity. Panels a and b in Figure 4 show the extracellular HA and 5HT concentrations when HA is released from the HA varicosity. In these simulations, both receptors on the 5HT varicosity are kept on. The release rate of the HA neuron is increased dramatically for 3.5 s (black bar). When the HA autoreceptor and 5HT receptor on the HA varicosity are both turned off, HA rises quite a lot (maroon). When the HA autoreceptor is on, but the 5HT receptor on the HA varicosity is kept off, there is not as big an increase in HA concentration (red). As discussed in the previous section, this is due to the autoreceptor inhibiting HA release as HA rises. Finally, when both the HA autoreceptor and the 5HT receptor on the HA varicosity are turned on, the perturbation in HA concentration is much smaller and the HA concentration comes back down below baseline quickly and oscillates (yellow). In this case, not only is the autoreceptor inhibiting HA release as HA increases, but as 5HT decreases (see panel b), the 5HT receptor on the HA varicosity withdraws stimulation of HA release. The corresponding concentrations of 5HT are shown in panel b. The most regulation of 5HT is seen when both the autoreceptor and the 5HT receptor on the HA varicosity are on (purple). This is because the corresponding HA concentrations do not rise as much (yellow). We explain the oscillations after the discussion of panels c and d.

Figure 4. Comodulation of 5HT and HA. Extracellular HA and 5HT concentrations after HA (panels a and b) or 5HT (panels c and d) stimulation. In panels a and b, the HA neuron is stimulated for 3.5 s (black bar). In these panels, both receptors on the 5HT varicosity are kept on. We show cases where both the HA autoreceptor and 5HT receptor on the HA varicosity are off (maroon and blue), the HA autoreceptor is on, but the 5HT receptor on the HA varicosity is not (red and green), and finally when both the HA autoreceptor and the 5HT receptor on the HA varicosity are on (yellow and purple). Similarly, in panels c and d the 5HT neuron is stimulated for 1.5 s (black bar). In these panels, both receptors on the HA varicosity are kept on. We show cases where both the 5HT autoreceptor and HA receptor on the 5HT varicosity are off (maroon and blue), the 5HT autoreceptor is on, but the HA receptor on the 5HT varicosity is not (red and green), and finally when both the 5HT autoreceptor and the HA receptor on the 5HT varicosity are on (yellow and purple). When all the receptors are on, each neurotransmitter adds an extra regulation of the other neurotransmitter.

Panels c and d in Figure 4 show the levels of HA and 5HT when the 5HT varicosity releases 5HT for 1.5 s (black bar). In these simulations, both receptors on the HA varicosity are kept on. Similar to our discussion above, the most regulation of 5HT and hence HA concentration is seen when both the 5HT autoreceptor and the HA receptor on the 5HT varicosity are on.

In both cases, oscillations arise when both receptors are on the respective varicosities. The reason for this stems from the comodulation between the two neurotransmitters. For example, in panels c and d, after stimulation of the 5HT neuron, the 5HT increases until 7 s and then stays above baseline until 9 s. This causes the HA to increase until 9 s. Then, since the HA stays above baseline until 11 s, the 5HT comes down below baseline causing the 5HT receptor on the HA varicosity to inhibit HA release, making HA go below baseline. However, as the HA goes below baseline, the HA receptor on the 5HT varicosity stimulates 5HT release causing 5HT to increase back to its baseline. This relationship continues for the rest of the 30 s period. It should be noted that these oscillations of HA are seen in our collaborator Parastoo Hashemi's lab as well; see Figure 8 in [9]. In conclusion, we have shown that the result of the comodulation between HA and 5HT is that each varicosity acts like an extra autoreceptor on the other varicosity, resulting in stronger regulation of the concentration of each neurotransmitter in the extracellular space.

4. Comodulation of histamine, dopamine, GABA and acetylcholine in the striatum

In the previous section, we showed that the comodulation of histamine and serotonin helps to stabilize both concentrations in the extracellular space. Conceptually, the ideas and results were fairly simple and straightforward. But, what are we to make of the situation depicted in Figure 5? Histamine (HA) neurons, dopamine (DA) neurons, Gamma-Aminobutyric Acid (GABA) neurons and acetylcholine (ACh) neurons project to the striatum, a region of the basal ganglia that controls motor movements (and other things). So, HA, DA, GABA and ACh are in the extracellular space in the striatum and each neurotransmitter affects the rate of release of other transmitters though numerous different receptors on the other neurons. These influences are inhibitory (red) or excitatory (green) as indicated in Figure 5; references are given in the figure legend. The challenge is to use mathematics to understand the dynamics of the whole system and its functions.

Figure 5. Comodulation in the Striatum. Varicosities are represented by the big rectangles and extracellular concentrations of neurotransmitters are represented by the small grey rectangles. HA modulates DA and GABA release via $H_3$ and $H_2$ receptors that inhibit release. HA modulates ACh release via $H_2$ and $H_1$ receptors that stimulate release and an $H_3$ receptor that inhibits release. DA modulates ACh release via a $D_1$ stimulatory receptor and $D_2$ inhibitory receptor and GABA release via a $D_1$ stimulatory receptor [33]. ACh modulates DA and GABA release via muscarinic receptors (mAChRs) that inhibit release and nicotinic receptors (nAChRs) that stimulate release [4]. GABA modulates ACh, HA and DA, release via inhibitory $GAB{A}_A$ receptors [33, 36]. In addition to having heteroreceptors, each varicosity has autoreceptors that inhibit release [8, 20, 26, 33].

For a number of reasons, this is challenging, indeed. For each receptor on a particular type of neuron, one would like to know the neurotransmitter concentration range over which the neurotransmitter affects release. One would like to know the maximal possible effects and, if possible, the kinetics. If one knows the density of that receptor, one could make a good mathematical model for the effect. In addition, total release will depend on the number of varicosities for each neurotransmitter. Unfortunately, for many receptors only partial information is known for rats or mice in particular tissues, and it can be problematic to simply assume that what is understood for rats and mice can be assumed for humans. Understanding the dynamics of this kind of complicated brain chemistry is important for understanding human behaviour in health and disease, but very challenging. In this section, we consider four different experimental situations in [32] and show how various subnetworks of the large network in Figure 5 give rise to the four sets of experimental results. This gives some understanding of the reasons for the behaviours seen in the experiments, but does not address the question of the dynamics of the whole network.

4.1. ACh release rate saturates as HA concentrations rise

HA has inhibitory and stimulatory receptors on ACh varicosities that affect the release of ACh [33]. In addition, there is a stimulatory ACh receptor on DA varicosities and an inhibitory DA receptor on acetylcholine varicosities [1]. We also take into account the DA autoreceptors present as in previous publications [8]. See the upper right diagram in Figure 6. Figure 6 shows a low (blue), medium (red) and high (green) level of HA administered from 20 to 30 min. The markers represent data redrawn from Figure 1 in [32]. Error bars are omitted but range from 0 to 0.4 units, and our model curves fall within them. Different levels of HA were modelled by increasing the input in the basal release function in the HA differential equation. Increased HA causes greater release of ACh, but this result saturates as higher levels of HA are given. This can be explained by the HA receptors present on the ACh varicosities. As HA levels rise, the receptors become saturated leading to very little changes in the ACh release rate. The role of the inhibitory HA and DA receptor is to slow down the rise of ACh release rates as HA increases. The DA autoreceptor causes the DA to come back down to steady state a bit faster so that the ACh release rates stay higher longer.

Figure 6. HA effect on ACh. ACh release rates are computed with a high (green), medium (red) and low (blue) level of histamine given from 20 to 30 min. The diagram in the upper right shows what receptors were included in the programme. The markers represent data redrawn from Figure 1 in [32]. Error bars are omitted but range from 0 to 0.4 units, and our model curves fall within them.

4.2. Decreased HA synthesis decreases ACh release rates

To analyze how extracellular HA concentrations affect ACh release, a HA synthesis inhibitor is administered, which reduces the amount of HA available for release. In this system, we only take into account the stimulatory HA receptor on ACh varicosities. See upper left diagram in Figure 7. To model the HA synthesis inhibitor α-fluoromethylhistidine, FMH, we decrease the basal release function in the HA differential equation from 20 to 30 min. See Section 5 for details of the model used. Figure 7 shows the ACh release rate when FMH is given from 20 to 30 min. The markers represent data redrawn from Figure 2 in [32]. Error bars on the markers range from 0.1 to 0.2 units, and our model curves fall within them. The ACh release rate decreases as HA synthesis decreases due to less stimulation on ACh release. This verifies the existence of the HA stimulatory receptor. Also, notice that the decrease in ACh release stabilizes in the last 20 min. Since the HA input stays at a decreased rate throughout the simulation, ACh is finding a new steady state.

Figure 7. Effect of an HA synthesis inhibitor. ACh release rates are computed when nothing is given (blue) or a HA synthesis inhibitor (red) α-fluoromethylhistidine, FMH, is given from 20 to 120 min. The diagram in the upper left shows what receptors were included in the programme. The markers represent data redrawn from Figure 2 in [32]. Error bars on the markers range from 0.1 to 0.2 units, and our model curves fall within them.

4.3. HA-induced ACh release rate decreases with DA receptor agonists/antagonists

ACh has $D_1$-like stimulatory DA receptors and $D_2$,$D_3$-like inhibitory DA receptors [33]. Quinpirole is a $D_2$,$D_3$ agonist and SKF 83566 is a $D_1$ antagonist. Thus quinpirole causes more inhibition of ACh release and SKF 83566 causes less stimulation of ACh release.

In Figure 8, we take into account these two DA receptors as well as the stimulatory HA receptor on ACh varicosities, the HA and DA autoreceptor, and the ACh stimulatory receptor on DA varicosities. See the diagram in the upper right. Three cases are observed: only HA is administered (blue), quinpirole, SKF 83566, and HA are administered (red) and only quinpirole and SKF 83566 are given (green). HA is given from 20 to 30 min and quinpirole and SKF 83566 are given 20 min before the simulation starts. Thus, for the red and green curves, the steady states for ACh, HA and DA are computed with quinpirole and SKF 83566 before simulation of the 70 min period.

To model HA administration, we increase the basal release in the HA differential equation from 20 to 120 min. To model quinpirole and SKF 83566 administration, the factor on the stimulatory DA receptor on the ACh varicosity is decreased and is increased on the inhibitory DA receptor on the ACh varicosity and the DA autoreceptor. Data is redrawn from Figure 3 in [32] in each case. The error bars are omitted due to space. However, our model curves fall within them. This figure shows that HA-induced ACh release rates decrease when quinpirole and SKF 83566 are given. The ACh release rate still rises as HA is given in the latter system, but it does not rise as much as when HA is given only. When histamine is given, the ACh release rate increases thus increasing the DA release rate. Since the ACh receptor on the DA varicosity is set to be stronger than the DA autoreceptor, DA concentrations increase. DA both inhibits and excites ACh release, but in the presence of quiniprole and SKF 83566, the excitation is inhibited and the inhibition becomes stronger. Thus, the ACh release rate increases (red), but not as much as the case without quinoprole and SKF 83566 (blue).

Figure 8. Effect of DA receptor agonists/antagonists. ACh release rates are computed when given only HA (blue), HA, quinpirole and SKF 83566 (red) and only quinpirole and SKF 83566 (green). Quinpirole is a $D_2$,$D_3$ DA receptor agonist (more inhibition of ACh) and SKF 83566 is a $D_1$ DA receptor antagonist (less stimulation of ACh). Quinpirole and SKF 83566 are given 20 min before timing for the 70 min period starts and HA is given from 20 to 30 min. The markers represent data redrawn from Figure 3 in [32]. Error bars are omitted but range from 0.1 to 0.3 units, and our model curves fall within them.

4.4. GABA receptor antagonist increases ACh release rate

An inhibitory GABA receptor is present on ACh varicosities. In Figure 9, a GABA receptor antagonist, bicuculline, is given from 20 to 70 min. We include the HA stimulatory receptor on ACh varicosities and the GABA inhibitory receptors on ACh and HA varicosities; see diagram in upper left. The HA stimulatory receptor was included so the effect of the GABA receptor antagonist is more prominent. The markers represent data redrawn from Figure 4 in [32]. Error bars are omitted but range from 0.1 to 0.4 units, and our model curves fall within them. To model bicuculine administration, the factors on the GABA receptors are decreased from 20 to 70 min. It is seen that as bicuculline is given the ACh release rate increases dramatically from 20 to 25 min, but starts to saturate in the last 40 min. This can be explained by ACh finding a new steady-state concentration to correspond to the lower inhibition via GABA receptors.

Figure 9. Effect of GABA receptor antagonist. ACh release rates are computed when a GABA receptor antagonist, bicuculline, is given from 20 to 70 min. The diagram shows the subnetwork modelled. The markers represent data redrawn from Figure 4 in [32]. Error bars are omitted but range from 0.1 to 0.4 units, and our model curves fall within them.

Figure 10. Schematic of the model. The figure shows the reactions in the model. The rectangular boxes indicate substrates and blue ellipses contain the acronyms of enzymes or transporters. Names of substrates are as follows: bHT, blood histidine; cHT, cytosolic histidine; HTpool, the histidine pool; cHA, cytosolic histamine; vHA, vesicular histamine; eHA, extracellular histamine; gHA, glial cell histamine; bTrp, blood tryptophan; Trp, cystolic tryptophan; Trp-pool, the tryptophan pool; 5-HTP, 5-Hydroxyindoleacetic acid; c5HT, cytosolic serotonin; v5HT, vesicular serotonin; e5HT, extracellular serotonin; g5HT, glial cell serotonin. Names of enzymes and transporters are as follows: Trpin, neutral amino acid transporter; TPH, tryptophan hydroxylase; AADC, aromatic amino acid decarboxylase; MAT, vesicular monoamine transporter; SERT, 5-HT reuptake transporter; MAO, monoamine oxidase; ALDH, aldehyde dehydrogenase; HTL, the histidine transporter; HTDC, histidine decarboxylase; HNMT, histamine methyltransferase; HAT, the putative HA transporter; H3, histamine autoreceptor (on HA varicosity) and heteroreceptor (on serotonin varicosity); R5HT, serotonin heteroreceptor; 5HT1b, serotonin autoreceptor.

5. Discussion

It has been known for a long time [17, 19] that volume transmission is an important phenomenon in the brain, but it is relatively unexplored in the mathematical biology community. In Section 2 we reviewed two homeostatic mechanisms in volume transmission. In ‘passive stabilization’ we explain why the concentration in the extracellular space of the projection region is relatively insensitive to cell death in the cell body region; an example is the concentration of DA in the striatum in the presence of cell death in the SNc in Parkinson's disease [5, 34]. Next, we showed that autoreceptors of varicosities provide end-product inhibition to synthesis, release and reuptake pathways, thus reducing the concentration changes in the projection region when there is a perturbation of firing rate.

In Section 3 we consider an example where two neurotransmitters in the extracellular space affect the rates of release of each other. HA inhibits 5HT release and 5HT stimulates HA release [7, 9, 15, 27, 29, 37]. We showed that this comodulation adds an extra regulation to each neurotransmitter ensuring that responses of concentrations to changes in firing rates are even smaller than the ones produced by autoreceptors. Essentially, each varicosity acts as an extra autoreceptor on the other varicosity, so the purpose of this comodulation is clear. Our computations were based on previously published models of HA and 5HT varicosities modified to include the comodulation effects and some other minor changes in parameters. The fact that each neurotransmitter acts like an extra autoreceptor on the release of the other neurotransmitter gives a clear purpose for this coregulation.

In Section 4, we confront a much harder situation where four neurotransmitters in the extracellular space, ACh, HA, DA and GABA all modulate each other in the striatum; see Figure 5. What is the purpose of this comodulation? To understand this requires investigation of a very complicated dynamical system and, unfortunately, only very partial information is available on the kinetics of the influences. So, as a beginning we consider four different experiments in Prast et al. [32, 33] and show that mathematical models of four different subnetworks give the behaviour seen in the experiments. These four (relatively simple) mathematical models are shown in Section 6.2.

The mathematical models used in Sections 3 and 4 have been investigated by machine computation, but mathematical analysis including steady-state analysis and bifurcation theory could shed light onto system behaviour and suggest new experiments. This is work we plan to do in the future.

Figure 5 was a diagram of local comodulation of neurotransmitters in the striatum, but it is likely that other similar diagrams exist in all major brain regions, perhaps with other neurotransmitters. The investigation of the dynamics of these diagrams is a major challenge for mathematicians and likely very important for understanding brain function in health and disease. For example, it is tempting to speculate that bipolar disorder and schizophrenia result from a bistable system with hysteresis. Whether this is true or not remains to be proven.

But there is a second major challenge for mathematicians. How does all of this volume transmission affect and regulate electrophysiological networks in which neurons make one-to-one connections with other neurons? For example, it is known that 5HT in the striatum binds to receptors on DA synapses on medium spiny neurons and increases the release of DA per action potential [11, 12]. Currently, in most electrophysiological models, all of this volume transmission is hidden in the parameters of the model. Thus, the interaction between the biochemical volume transmission system and the electrophysiological networks raises open, compelling questions. For examples of work that addresses these issues in particular systems, see [13, 31].

Although we have written as though the HA neurons, DA neurons, GABA neurons and ACh neurons are distinct types of neurons, it is now widely accepted that many neurons corelease multiple neurotransmitters [23, 39]. In fact, there is evidence that the dopaminergic neurons of the SNc corelease GABA, with consequences for the direct pathway of the basal ganglia and for behaviour [21]. Our model is actually agnostic concerning the neuronal origin of some of the neurotransmitters considered (e.g. dopamine and GABA), while we do assume that HA and 5HT are released from different cells. Neurotransmitter corelease raises additional questions that we will address in ongoing work.

6. Models

In this section, we describe the mathematical models used to produce the results found in this paper. For Figures 2– 4 we combined our previous histamine and serotonin models, simplified the existing receptors, and added a serotonin receptor on the histamine varicosity [7, 9]. Our combined programme is shown schematically in Figure 10. We discuss other small changes to the integrated model in Section 6.1. For Section 4, we model different subnetworks of Figure 5 for each experimental situation. Section 6.2 provides the models used for the figures in Section 4 (Figures 6– 9) including the differential equations and the kinetics of receptors assumed.

6.1. Serotonin and histamine comodulation

When combining our previous serotonin and histamine programmes we simplified the models a little. For both varicosities, we simplified the receptor and autoreceptor models by allowing the strength of the receptor effects to depend simply on the current concentrations in the extracellular space. These receptor functions can be found in Table 1. We also removed dihydrobiopterin ($bh2$), tetrahydrobiopterin ($bh4$) and 5-hydroxyindoleacetic acid (hiaa) because they play no role in the current simulations. Eliminating these variables from the differential equations in [7, 9] and adding the serotonin receptor on the histamine varicosity gives the differential equations used in this paper.

Table 1. Functions for Receptors (μM,μM/hr,/hr).

Receptor	Function
HA autoreceptor	$=1-\left(0.5\right)\cdot {\displaystyle {\displaystyle \frac{eHA-eH{A}_{ss}}{1+|eHA-eH{A}_{ss}|}}}$
5HT receptor on HA var.	${\displaystyle =10+\left(15\right)\cdot {|2\cdot (e5HT-e5H{T}_{ss}\left)\right|}^{{\displaystyle \frac{1}{3}}}}$
5HT autoreceptor	$=1.89-\left(20\right)\cdot {\displaystyle {\displaystyle \frac{e5HT-e5H{T}_{ss}}{1+|e5HT-e5H{T}_{ss}|}}}$ (Figure 3)
	$=1.89-\left(4\right)\cdot {\displaystyle {\displaystyle \frac{e5HT-e5H{T}_{ss}}{1+|e5HT-e5H{T}_{ss}|}}}$ (Figure 4)
HA receptor on 5HT var.	$=1-{\displaystyle {\displaystyle \frac{eHA-eH{A}_{ss}}{1+|eHA-eH{A}_{ss}|}}}$

Figure 10 is the schematic diagram that was used for our integrated programme. For each varicosity, the model includes synthesis, storage, release, reuptake and removal of the neurotransmitter. For the HA terminal, histidine in the blood (bHT) is transported into the varicosity by the histidine transporter (HTL) where it becomes cytosolic histidine (cHT) or goes into the histidine pool (HTpool). Most of the histidine that enters the cell is used for other processes than making HA and that is what the HTpool represents. cHT is converted to cytosolic histamine, cHA, by the enzyme histidine decarboxylase, HTDC. Some cHA is catabolized by the enzyme histamine methyltransferase, HNMT, some is transported into the vesicles by the monoamine transporter, MAT, and becomes vesicular HA, vHA, and some leaks out of the cytosol into the extracellular space (indicated by the dashed line). vHA is released into the extracellular space, at a rate proportional to neuronal firing, where it becomes extracellular HA, eHA. In the extracellular space, eHA has several fates. It can be transported back into the cytosol by a putative HA transporter, HAT. Although HAT has not yet been identified, it must exist since unpublished experimental data (personal communication [22]) shows that it is hard to deplete vesicle stores and this suggests strongly that HA must be transported back into the varicosity from the cytosol. It can diffuse away (removal). It can be transported into glial cells where it becomes glial HA, gHA, which then leaks out or is catabolized by HNMT [18, 24, 40]. Finally, eHA is regulated by the autoreceptor H3 and the serotonin R5HT receptor. Similarly, for the 5HT terminal, the model includes uptake of tryptophan (bTrp) across the blood-brain barrier and transport into terminals; synthesis of 5HT by tryptophan hydroxylase (TPH) and aromatic amino acid decarboxylase (AADC); transport of 5HT into a vesicular compartment by the monoamine transporter (MAT); release of 5HT into the extracellular space depending on the firing rate; reuptake via the SERTs; regulation by the autoreceptor 5HT1b and histamine H3 receptor; and diffusion away from the terminal defined as removal.

Some changes were made to the coefficients and kinetic parameters for the enzymes from [7, 9]. We save the receptor functions in Table 1. For the histamine varicosity, the coefficients that were changed and their new values are: leakage coefficient of histamine from the cytosol to the extracellular space is 12; histamine release per action potential is 0.82; removal of extracellular histamine is 4; transfer of the cytosolic histamine pool to cytosolic histamine is 0.1. The enzyme rates that were changed and their new parameters are: HTDC, $V_{max}:117$; MAT, $V_{max}:6300$; HAT, $V_{max}:1172.34$; HAT (glial), $V_{max}:5000$, $k_m:20$. For the serotonin varicosity, all coefficients remain the same. The enzyme rates that were changed and their new parameters are: TPH, $V_{max}$: 27.8; AADC, $V_{max}$: 1200; CATAB, $V_{max}$: 8000.

For passive stabilization (Figure 2), we decoupled the serotonin and histamine varicosities by setting the histamine H3 receptor on the serotonin varicosity to have zero effect on serotonin, so we could analyze the serotonin varicosity separately. To implement the effect of serotonin neurons dying, we multiplied all of the serotonin equations by a factor f to represent the fraction of cells alive except the differential equation for extracellular serotonin. For this equation, we multiplied the release and reuptake of serotonin by f and divided the removal term by f, because removal increases as there are fewer and fewer varicosities with SERTs available. We, also, decreased the coefficient of the extracellular serotonin removal term to 1. Figures 3 and 4 used the coupled serotonin and histamine ODE system. Subscript ss represents concentration at steady state.

6.2. Comodulation of HA, DA, GABA and ACh in the striatum

Below we give the model used for each experimental situation in Figures 6– 9. In these figures, we modelled the increase/decrease of release of ACh. The actual drug concentrations or doses were not used. $basalre{l}_X$ is the basal release rate and $d_X$ is the catabolism coefficient of neurotransmitter X. All basal release rates are $=10\mu M/hr$ except for ACh in Figure 8 which is 1. Each figure indicates the subnetwork that we are assuming. For each variable X, the rate of change is $basalre{l}_X$ times the receptor functions for other variables that influence the release minus the catabolism rate. The receptor functions for each neurotransmitter are defined as $R{X}_{Z+/-}$, the influence of the Z stimulatory/inhibitory receptor on the varicosity of X. They multiply the release rate by 1 plus/minus a Michaelis Menten function depending on whether they are stimulatory or inhibitory, respectively. The differential equations for the four cases follow. Table 2 provides the receptor functions, basal release rates and catabolism rates for each experimental situation.

Table 2. Functions and Coefficients (μM,μM/hr,/hr).

Velocity	Parameter	Model Value
Figure 6
$RD{A}_{\text{ ACh+}}$	ACh + rec. on DA	if ACh $<eq.$: $=1+{\displaystyle \frac{V_{max}(ACh-AC{h}_{ss})}{k_m+(ACh-AC{h}_{ss})}}$
		if ACh $>eq.$: $=1-{\displaystyle \frac{V_{max}(ACh-AC{h}_{ss})}{k_m+(ACh-AC{h}_{ss})}}$
	$k_m$	1
	$V_{max}$	15
$RD{A}_{\text{ DA-}}$	DA Autoreceptor	if DA $<eq.$: $=1+{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
		if DA $>eq.$: $=1-{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
	$k_m$	0.0001
	$V_{max}$	0.9
$RAC{h}_{\text{ HA+}}$	HA + rec. on ACh	if HA $<eq.$: $=1-{\displaystyle \frac{V_{max}(HA-H{A}_{ss})}{k_m+(HA-H{A}_{ss})}}$
		if HA $>eq.$: $=1+{\displaystyle \frac{V_{max}(HA-H{A}_{ss})}{k_m+(HA-H{A}_{ss})}}$
	$k_m$	50
	$V_{max}$	10
$RAC{h}_{\text{ HA-}}$	HA – rec. on ACh	if HA $<eq.$: $=1+{\displaystyle \frac{V_{max}(HA-H{A}_{ss})}{k_m+(HA-H{A}_{ss})}}$
		if HA $>eq.$: $=1-{\displaystyle \frac{V_{max}(HA-H{A}_{ss})}{k_m+(HA-H{A}_{ss})}}$
	$k_m$	10
	$V_{max}$	0.6
$RAC{h}_{\text{DA-}}$	DA – rec. on ACh	if DA $<eq.$: $=1+{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
		if DA $>eq.$: $=1-{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
	$k_m$	0.5
	$V_{max}$	0.8
	$d_1$	2
	$d_2$	2
	$d_3$	2
	$d_4$	2
Figure 7
$RAC{h}_{\text{HA+}}$	HA + rec. on ACh	$=1+{\displaystyle \frac{V_{max}HA}{k_m+HA}}$
	$k_m$	12
	$V_{max}$	3
Catabolism Rates
	$d_1$	2
	$d_2$	2
	$d_3$	2
	$d_4$	2
Velocity	Parameter	Model Value
Figure 8
$RH{A}_{\text{ HA-}}$	HA autoreceptor	if HA $<eq.$: $=1+{\displaystyle \frac{V_{max}(HA-H{A}_{ss})}{k_m+(HA-H{A}_{ss})}}$
		if HA $>eq.$: $=1-{\displaystyle \frac{V_{max}(HA-H{A}_{ss})}{k_m+(HA-H{A}_{ss})}}$
	$k_m$	7
	$V_{max}$	0.8
$RD{A}_{\text{ ACh+}}$	ACh rec. on DA	if ACh $<eq.$: $=1-{\displaystyle \frac{V_{max}(ACh-AC{h}_{ss})}{k_m+(ACh-AC{h}_{ss})}}$
		if ACh $>eq.$: $=1+{\displaystyle \frac{V_{max}(ACh-AC{h}_{ss})}{k_m+(ACh-AC{h}_{ss})}}$
	$k_m$	0.7
	$V_{max}$	10
$RD{A}_{\text{ DA-}}$	DA Autoreceptor	if DA $<eq.$: $=1+{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
		if DA $>eq.$: $=1-{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
	$k_m$	0.0001
	$V_{max}$	0.7 (.999 for Quin. and SKF 83566)
$RAC{h}_{\text{ HA+}}$	HA + rec. on ACh	if HA $<eq.$: $=1-{\displaystyle \frac{V_{max}(HA-H{A}_{ss})}{k_m+(HA-H{A}_{ss})}}$
		if HA $>eq.$: $=1+{\displaystyle \frac{V_{max}(HA-H{A}_{ss})}{k_m+(HA-H{A}_{ss})}}$
	$k_m$	100
	$V_{max}$	3.5
$RAC{h}_{\text{ DA+}}$	DA + rec. on ACh	if DA $<eq.$: $=1-{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
		if DA $>eq.$: $=1+{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
	$k_m$	0.0001
	$V_{max}$	1 (0.01 for Quin. and SKF 83566)
$RAC{h}_{\text{ DA-}}$	DA – rec. on ACh	if DA $<eq.$: $=1+{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
		if DA $>eq.$: $=1-{\displaystyle \frac{V_{max}(DA-D{A}_{ss})}{k_m+(DA-D{A}_{ss})}}$
	$k_m$	1.5
	$V_{max}$	0.8 (.9999 for Quin. and SKF 83566)
Catabolism Rates
	$d_1$	10
	$d_2$	2
	$d_3$	2
	$d_4$	2
Figure 9
$RH{A}_{\text{ GABA-}}$	GABA – rec. on HA	if t<20 min: $=1-{\displaystyle \frac{V_{max}GABA}{k_m+GABA}}$
		if t>20 min: increases to 0.99 linearly
	$k_m$	3
	$V_{max}$	0.9
$RAC{h}_{\text{ HA+}}$	HA + rec. on ACh	$=1+{\displaystyle \frac{V_{max}HA}{k_m+HA}}$
	$k_m$	7
	$V_{max}$	3
$RAC{h}_{\text{ GABA-}}$	GABA – rec. on ACh	if t<20 min: $=1-{\displaystyle \frac{V_{max}GABA}{k_m+GABA}}$
		if t>20 min: increases to 0.99 linearly
	$k_m$	3
	$V_{max}$	0.8
Catabolism Rates
	$d_1$	2
	$d_2$	2
	$d_3$	2
	$d_4$	2

Figure 6 Differential Equations $\begin{array}{rl}\frac{\mathrm{d}\left[HA\right]}{\mathrm{d}t}& =basalre{l}_{\text{HA}}\left(t\right)-d_1\cdot HA\\ \frac{\mathrm{d}\left[GABA\right]}{\mathrm{d}t}& =basalre{l}_{\text{GABA}}\left(t\right)-d_2\cdot GABA\\ \frac{\mathrm{d}\left[DA\right]}{\mathrm{d}t}& =basalre{l}_{\text{DA}}\left(t\right)\cdot RD{A}_{\text{ACh+}}\left(ACh\right)\cdot RD{A}_{\text{DA-}}\left(DA\right)-d_3\cdot DA\\ \frac{\mathrm{d}\left[ACh\right]}{\mathrm{d}t}& =basalre{l}_{\text{ACh}}\left(t\right)\cdot RAC{h}_{\text{HA+}}\left(HA\right)\cdot RAC{h}_{\text{HA-}}\left(HA\right)\\ & \times RAC{h}_{\text{DA-}}\left(DA\right)-d_4\cdot ACh\end{array}$Figure 7 Differential Equations $\begin{array}{rl}\frac{\mathrm{d}\left[HA\right]}{\mathrm{d}t}& =basalre{l}_{\text{HA}}\left(t\right)-d_1\cdot HA\\ \frac{\mathrm{d}\left[GABA\right]}{\mathrm{d}t}& =basalre{l}_{\text{GABA}}\left(t\right)-d_2\cdot GABA\\ \frac{\mathrm{d}\left[DA\right]}{\mathrm{d}t}& =basalre{l}_{\text{DA}}\left(t\right)\cdot RD{A}_{\text{DA-}}\left(DA\right)-d_3\cdot DA\\ \frac{\mathrm{d}\left[ACh\right]}{\mathrm{d}t}& =basalre{l}_{\text{ACh}}\left(t\right)\cdot RAC{h}_{\text{HA+}}\left(HA\right)-d_4\cdot ACh\end{array}$Figure 8 Differential Equations $\begin{array}{rl}\frac{\mathrm{d}\left[HA\right]}{\mathrm{d}t}& =basalre{l}_{\text{HA}}\left(t\right)\cdot RH{A}_{\text{HA-}}\left(HA\right)-d_1\cdot HA\\ \frac{\mathrm{d}\left[GABA\right]}{\mathrm{d}t}& =basalre{l}_{\text{GABA}}\left(t\right)-d_2\cdot GABA\\ \frac{\mathrm{d}\left[DA\right]}{\mathrm{d}t}& =basalre{l}_{\text{DA}}\left(t\right)\cdot RD{A}_{\text{ACh+}}\left(ACh\right)\cdot RD{A}_{\text{DA-}}\left(DA\right)-d_3\cdot DA\\ \frac{\mathrm{d}\left[ACh\right]}{\mathrm{d}t}& =basalre{l}_{\text{ACh}}\left(t\right)\cdot RAC{h}_{\text{HA+}}\left(HA\right)\\ & \times RAC{h}_{\text{DA+}}\left(DA\right)\cdot RAC{h}_{\text{DA-}}\left(DA\right)-d_4\cdot ACh\end{array}$Figure 9 Differential Equations $\begin{array}{rl}\frac{\mathrm{d}\left[HA\right]}{\mathrm{d}t}& =basalre{l}_{\text{HA}}\left(t\right)\cdot RH{A}_{\text{GABA-}}\left(GABA\right)-d_1\cdot HA\\ \frac{\mathrm{d}\left[GABA\right]}{\mathrm{d}t}& =basalre{l}_{\text{GABA}}\left(t\right)-d_2\cdot GABA\\ \frac{\mathrm{d}\left[DA\right]}{\mathrm{d}t}& =basalre{l}_{\text{DA}}\left(t\right)-d_3\cdot DA\\ \frac{\mathrm{d}\left[ACh\right]}{\mathrm{d}t}& =basalre{l}_{\text{ACh}}\left(t\right)\cdot RAC{h}_{\text{HA+}}\left(HA\right)\cdot RAC{h}_{\text{GABA-}}\left(GABA\right)-d_4\cdot ACh\end{array}$

Acknowledgments

The authors thank Parry Hashemi, Imperial College, London, who brought the idea of comodulation to our attention and provided references.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Allison Cruikshank (AC) was supported by National Science Foundation DMS-2038056, a research and training grant (PI Alexander Kiselev, Duke University).