Comodulation of dopamine and serotonin on prefrontal cortical rhythms: a theoretical study

The prefrontal cortex (PFC) is implicated to play an important role in cognitive control. Abnormal PFC activities and rhythms have been observed in some neurological and neuropsychiatric disorders, and evidences suggest influences from the neuromodulators dopamine (DA) and serotonin (5-HT). Despite the high level of interest in these brain systems, the combined effects of DA and 5-HT modulation on PFC dynamics remain unknown. In this work, we build a mathematical model that incorporates available experimental findings to systematically study the comodulation of DA and 5-HT on the network behavior, focusing on beta and gamma band oscillations. Single neuronal model shows pyramidal cells with 5-HT1A and 2A receptors can be non-monotonically modulated by 5-HT. Two-population excitatory-inhibitory type network consisting of pyramidal cells with D1 receptors can provide rich repertoires of oscillatory behavior. In particular, 5-HT and DA can modulate the amplitude and frequency of the oscillations, which can emerge or cease, depending on receptor types. Certain receptor combinations are conducive for the robustness of the oscillatory regime, or the existence of multiple discrete oscillatory regimes. In a multi-population heterogeneous model that takes into account possible combination of receptors, we demonstrate that robust network oscillations require high DA concentration. We also show that selective D1 receptor antagonists (agonists) tend to suppress (enhance) network oscillations, increase the frequency from beta toward gamma band, while selective 5-HT1A antagonists (agonists) act in opposite ways. Selective D2 or 5-HT2A receptor antagonists (agonists) can lead to decrease (increase) in oscillation amplitude, but only 5-HT2A antagonists (agonists) can increase (decrease) the frequency. These results are comparable to some pharmacological effects. Our work illustrates the complex mechanisms of DA and 5-HT when operating simultaneously through multiple receptors.


INTRODUCTION
The prefrontal cortex (PFC) plays an essential role in many higher brain functions such as goal-directed behavior, action planning, learning, attention, mnemonic processes, inhibitory control, and task switching (Miller, 2000;Fuster, 2001;Miller and Cohen, 2001;Andrade, 2011b). Neural oscillations in the PFC are suggested to be important for communication within the PFC and with other brain regions, and are suggested to regulate such higher cognitive functions (Wang, 2010;Benchenane et al., 2011).
At the neuronal network level, it has been found that DA injected in the PFC of anesthetized rats enhances hippocampalprefrontal coherence in the theta band oscillation (Benchenane et al., 2010), which could be due to DA modulating the GABAergic inhibition (Tierney et al., 2008). Blocking D1 receptors has been known to increase alpha and beta band oscillations more in local field potentials for novel than familiar associations (Puig and Miller, 2012). Increasing extracellular DA with genetic polymorphism of dopamine transporter (DAT1) in humans can enhance evoked gamma response to stimulus (Demiralp et al., 2007) 5-HT can also increase the frequency and amplitude of slow waves by promoting the UP states in PFC via activation of 5-HT2A receptors, suggesting an excitatory effect in in vivo condition (Puig et al., 2010). 5-HT2A/2C receptor agonist/antagonist has also been found to synchronize/desynchronize frontal cortical oscillations in anesthetized rats (Budzinska, 2009).
Although there have been extensive investigations on the modulation of DA and 5-HT on the PFC, little is known about their comodulation effects on the PFC network dynamics and their potential applications in drug treatments (Diaz-Mataix et al., 2005;Di Pietro and Seamans, 2007;Artigas, 2010). In fact, many of the DA and 5-HT induced intracellular signaling pathways overlap (Amargos-Bosch et al., 2004;Santana et al., 2004;Di Pietro and Seamans, 2007;Esposito et al., 2008;Santana et al., 2009), suggesting that DA and 5-HT may cooperatively modulate PFC activity. One notable study has found that coadministration of 5-HT2A antagonist with a D2 antagonist in PFC significantly increase DA release which is greater than that induced by either antagonist alone (Westerink et al., 2001). A recent research has found that co-application of DA and 5-HT can increase the evoked excitability of certain PFC pyramidal cells (the gain of the neuronal input-output response) more than when either was applied alone, while the activities of other pyramidal cells get more suppressed (Di Pietro and Seamans, 2011). Furthermore, the same study also shows that prior DA or 5-HT application can potentiate the subsequent effect of the other.
In this work, we integrate the essential available experimental findings into a biologically motivated computational model to provide insights into the possible PFC dynamics caused by the comodulation of DA and 5-HT. The focus will be on tonic DA and 5-HT modulations, and their effects on higher frequency band oscillations.

MATERIALS AND METHODS
The computational models in this work will implement DA and 5-HT comodulation at the neuronal and synaptic levels.
In the following, we shall discuss about the various neuronal constituents of the PFC modulated by DA and 5-HT.
The instantaneous population firing rate for each neuronal subgroup follows the established dynamics (Wilson and Cowan, 1972;Dayan and Abbott, 2001;Murphy and Miller, 2009): where i is the index for the subgroup of neurons. For pyramidal cells, i = 1 to 4 denote the subpopulation expressing D1+5HT1A, D2+5HT1A, D1+5HT1A+5-HT2A, and D2+5HT1A+5-HT2A, respectively. For interneurons, i = 5 to 8 denote the subgroup neurons expressing D1+5HT1A, D1+5HT2A, D2+5HT1A, and D2+5HT2A , respectively. τ i is the neuronal membrane time constant, set at 10 ms for pyramidal cells and 15 ms for inhibitory neurons (McCormick et al., 1985). I i, syn is the total synaptic currents to the i-th subgroup. The activation function f i (I i, syn ) follows the established form (Eckhoff et al., 2011): where r max denotes the saturated firing rate, r max is 80 Hz for pyramidal cells, and 120 Hz for inhibitory neurons. C is the  Gaspar et al., 1995;Amargos-Bosch et al., 2004;Santana et al., 2004Santana et al., , 2009Andrade, 2011a;Puig, 2011. gain and set as 300 Hz/nA for pyramidal cells and 500 Hz/nA for inhibitory interneurons. I L is associated with the membrane leakage current and its value is set at 150 Hz for pyramidal cells and 180 Hz for inhibitory interneurons. The curvature of the activation function g is 0.2 Hz −1 .

MODULATION OF NEURONAL EXCITABILITY
In DA modulation, although D1 and D2 receptors can act on different signaling transduction pathways, they can effectively attain opposite effects (e.g., D1 activate protein kinase A while D2 receptors inactivate it) Beaulieu and Gainetdinov, 2011;Tritsch and Sabatini, 2012). Therefore, we model the modulation of D1 (D2) on neuronal activity by increasing (decreasing) the gain factor C and decreasing (increasing) the leakage factor I L of the input-output function (Thurley et al., 2008). With regard to 5-HT modulation, experiments have shown that 5-HT1A can hyperpolarize the neurons through the activation of G protein-gated inwardly rectifying K + channels, while 5-HT2A activation can induce slow membrane depolarization and inhibition of the slow after-hyperpolarization which increases membrane excitability (Andrade et al., 1986;Beique et al., 2004;Andrade, 2011a). 5-HT1A receptors are often localized on the axon initial segment and soma of pyramidal neurons where they act to suppress action potential generation, while 5-HT2A receptors are abundant in apical dendrites where they can amplify the synaptic current (Amargos-Bosch et al., 2004;Santana et al., 2004). 5-HT2A receptors are also found to increase the gain of the input-output relationship of pyramidal neurons (Zhang and Arsenault, 2005). Therefore, the modulation of 5-HT on neuronal activity can be modeled by an increase in the leakage factor I L (for 5-HT1A) and increase of the gain factor C (for 5-HT2A).
The activation of D1, D2, 5-HT1A and 5-HT2A receptors are concentration dependent Hurley, 2006;Solt et al., 2007). For simplicity, we apply the sigmoid function to describe the concentration dependent modulation of DA and 5-HT on the gain and leak factors (Tables 2 and 3), similar to previous work (Fellous and Linster, 1988;Scheler, 2004). In these formulae, [DA] 1 denotes the half maximal effective concentration (EC 50 ) of DA for D1 receptor, and [DA] 2 is the EC 50 of DA for D2 receptor. Similarly,  1 and  2 are the EC 50 for 5-HT1A and 5-HT2A receptors, respectively. In brief, we depict the modulation of DA and 5-HT on neuronal activity by changing the activation function through multiplying C with a gain factor shown in Table 2, and I L with the leakage factor shown in Table 3.
The concentrations of DA and 5-HT in Tables 2 and 3 can be inferred from experiments such as those using microdialysis and voltammetry techniques. It is shown that the basal extracellular DA and 5-HT concentrations in the PFC is about ∼0.2-2.5 nM/L in resting condition, and can increase by as much as 10-200% when performing behavioral tasks (Adell et al., 1991;Watanabe et al., 1997;Lena et al., 2005;Winstanley et al., 2006;Rogoz and Golembiowska, 2010;Seeman, 2010;Staiti et al., 2011;van Dijk et al., 2012).

Neuron
Receptors Gain factor

See text for parameter values.
D1-and D2-like receptors can have a high affinity state with a binding constant around the nM/L level or a low affinity state with a binding constant around the μM/L level (Richfield et al., 1989). Since we are studying only tonic concentration levels, we shall only focus on the high affinity receptors, by assuming that low affinity ones are activated more in the phasic or evoked mode, e.g., during behavioral tasks. In particular, we assume that the high affinity D1-and D2-like receptors are sensitive within a range of 0-50 nM/L with different EC 50 values. We chose [DA] 1 = 4 nM/L and [DA] 2 = 8 nM/L (Koshkina, 2006), suggesting lower [DA] activates D1 receptor only while higher [DA] activates both D1 and D2 receptors, similar to the observed activation order of these receptors depending on DA concentration .
Similarly, 5-HT1A and 5-HT2A receptors can also operate in high-affinity and low-affinity states (Glennon et al., 1998;Watson et al., 2000). In this work, we assume that 5-HT1A and 5-HT2A receptors operate at high-affinity state since tonic 5-HT concentration at the nM/L level is far below the affinity of 5-HT for the low agonist affinity state (Watson et al., 2000), and we vary the 5-HT concentration within the range 0-5 nM/L. The affinity of 5-HT for 5-HT1 is higher than most of other subtype of

SYNAPTIC CURRENTS
The synaptic currents are mediated by AMPA, NMDA, and GABA receptors. We adopt an all-to-all connectivity among neuronal subgroups, thus, the synaptic currents to one neuron in the i-th subgroup can be approximated by summing all presynaptic neurons (N j ) and normalized by the neurons in the i-th subgroup (N i ): where the term τ A G A, ext, i r ext describes the constant part of the AMPA mediated background Poisson input with rate r ext (2.4 KHz) (Wong and Wang, 2006). Unlike (Wong and Wang, 2006;Eckhoff et al., 2011), we do not include synaptic noise into the model as we find that noise does not significantly affect our results (not shown). G X, ij in Equation (3) denotes the averaged coefficient or strength of the synaptic currents mediated by receptors X (A for AMPA, N for NMDA, and G for GABA) from neuron j to i. Their values are constrained by the experimentally observed neural circuit oscillation frequencies and are assigned as: G A, PP = 4.42 nA, G A, IP = 4.21 nA, G N, PP = 0.10 nA, G N, IP = 0.83 nA, G G, PI = 2.275 nA, G G, II = 1.75 nA, G A, ext, P = 0.0929 nA, and G A, ext, I = 0.0716 nA. S X, j is the averaged synaptic gating variable mediated by receptors X expressing on neuron in the j-th subgroup and follows the established dynamical forms (Brunel and Wang, 2001;Wong and Wang, 2006;Eckhoff et al., 2011): where τ A = 2 ms, τ N = 100 ms, and τ G = 10 ms are the decay time constants for AMPA, NMDA, and GABA receptors, respectively. The fraction N j /N i can be approximated according to the experimental observations on the DA and 5-HT receptors distribution (Pazos and Palacios, 1985;Amargos-Bosch et al., 2004;Beique et al., 2004;Santana et al., 2004;de Almeida and Mengod, 2007;Santana et al., 2009). Based on experimental observations (see Table 1), we assume that the ratio for pyramidal cells expressing D1 to D2 receptors is approximately 25:15%, and the ratio for interneurons expressing D1 to D2 receptors is approximately 30:10%. We also assume that the ratio for pyramidal cells or interneurons solely expressing 5-HT1A receptors to that those expressing 5-HT1A and 5-HT2A receptors is approximately 50:50%. Moreover, pyramidal cells expressing D2 receptors are often found in apposition with GABAergic cells not expressing D2 receptors, leading to the synaptic currents from inhibitory neurons expressing D2 receptors smaller than those from inhibitory neurons expressing D1 receptors (de Almeida and Mengod, 2007). So we specify the fraction from D2 expressing interneurons to D2 expressing pyramidal cells as a fifth of that from D1 expressing interneurons. Table 4 lists the above-mentioned fractions, where the fraction from the j-th to i-th subgroup is the value at the i-th row and j-th column. When simulating the two-population model, we ignore other subgroups by setting the irrelevant connections to zero.

MODULATION OF SYNAPTIC TRANSMISSION
The synaptic current coefficients or strengths, G X, ji , can be modulated by DA and 5-HT through the activation of D1, D2, 5-HT1A, and 5-HT2A receptors. Studies have demonstrated that D1-like receptors can enhance AMPA, NMDA, and GABA mediated synaptic currents or receptor expression, while D2-like receptors decrease them (Seamans et al., 2001b;Gorelova et al., 2002;Thurley et al., 2008). Similar to the modulation of DA on synaptic currents, 5-HT also bidirectionally modulates the synaptic currents through the activation of different receptors: 5-HT1A receptors reduce AMPA and NMDA mediated currents (Cai et al., 2002;Zhong et al., 2008), while 5-HT2A receptors increase them. Activation of 5-HT2A receptors on pyramidal cells can reduce GABA mediated currents (Feng et al., 2001) while 5-HT1A receptors may suppress the presynaptic GABAergic release in interneurons (Yan, 2002). We adopt a scalar, sigmoid function factor for the modulation of synaptic current coefficients (Fellous and Linster, 1988;Scheler, 2004). In principle, as all 4 considered receptors can modulate the three (AMPA, NMDA and GABA mediated) synaptic currents, we have 12 possible modulation factors ( Table 5). The parameters λ ij and κ ij depict the modulatory effects of the i-th receptor on the j-th synaptic type. We set the amplitude λ ij = 0.2 and the curvature κ ij = 1/nM for i = 1, 2 and κ ij = 4/nM for i = 3, 4.
Usually, the synaptic modulation factors are determined by the specific types of receptors expressing on the postsynaptic neuron. 1 + e −κ 43 ( −  2 ) . Figure 2A shows that it is more effective to have [DA] and  to increase/decrease in the same direction. One exception is the modulation of 5-HT1A on GABA-mediated synaptic currents. As mentioned above, 5-HT1A receptors can modulate the GABA-mediated currents via a presynaptic mechanism by reducing the GABAergic release (Yan, 2002), which effectively reduces the GABA-mediated currents. Thus, we assume that only presynaptic expressing 5-HT1A receptors can modulate GABAmediated currents by the factor 1 −  Figure 2B shows that having DA and 5-HT to increase/decrease in opposite directions provides a more effective modulation effect on this particular synapse.

RESULTS
In the following, we shall first investigate the effects of DA and 5-HT modulation on individual PFC neurons, and then followed by various coupled excitatory-inhibitory PFC circuits. After that, we investigate how a more realistic heterogeneous multi-population network model is modulated by DA and 5-HT concentration levels and selective receptor agonists/antagonists.

NON-MONOTONIC MODULATION OF 5-HT ON PYRAMIDAL CELL COEXPRESSING 5-HT1A AND 5-HT2A RECEPTORS
As previously mentioned, our model assumes that D1 and D2 receptors are expressed on distinct neuronal populations. Thus, it is expected that the modulation of DA on single neuronal activity should monotonically depend on the extracellular DA concentration. Similarly, for neurons solely expressing 5-HT1A or 5-HT2A receptors, the neuronal activity will also monotonically depend on extracellular 5-HT concentration. However, for pyramidal cells coexpressing 5-HT1A and 5-HT2A receptors, the modulation of 5-HT is not monotonic. The steady firing rate for such a pyramidal cell is r =  2 ) C P I syn − 1 +¯δ 1 1+e −ᾱ 1 ( −  1 ) I L . For any constant I syn , the combined modulation effect is found to be non-monotonic (Figure 3). Lower 5-HT concentration initially decreases the neuronal firing rate because of the activation of the high affinity of 5-HT1A receptors. But higher concentration of 5-HT subsequently increases the neuronal firing rate due to activating the lower affinity 5-HT2A receptors.

MODULATION ON TWO-POPULATION EXCITATORY-INHIBITORY NEURONAL NETWORKS
The simplest "canonical" cortical column and its oscillatory behavior can be modeled by an excitatory neuronal population mutually coupled to an inhibitory neuronal population (Wilson and Cowan, 1972). In our model, for every pair of coupled excitatory and inhibitory neuronal populations considered, we set the parameters (i.e., the fractions in Table 4) describing the other (six) populations to be zero. Based on the expression of the various receptors (Tables 2 or 3), there are 4 × 4 = 16 combinations of excitatory and inhibitory neurons in the two-population network model. In general, we find that if the network consists of pyramidal cells which express D2 receptors (Pyr2 or Pyr4), the network cannot attain oscillatory behavior over the ranges of [DA] and  explored. However, if the network includes D1-expressing pyramidal cells (Pyr1 or Pyr3), a rich repertoire of dynamical behavior can be produced with varying [DA] and , as described below. Figure 4 shows the results for the Pyr1-type (D1+5-HT1A) pyramidal cells coupled to Int1-type (D1+5HT1A) interneurons. Figure 4A shows an example of the firing rate time course of Pyr1 neurons for [DA] = 7 nM. Higher  level results in lower oscillation amplitude but faster frequency. The firing rate time courses for the inhibitory Int1 populations look similar (not shown). At intermediate [DA], the oscillation frequency decreases with increasing [DA] while remaining within the high beta range ( Figure 4B). Figure 5A1 (green lines) summarizes the oscillation amplitudes over a range of [DA] levels. The top green lines denote the maximum (top) and minimum (bottom) firing rates during oscillation. The red lines represent collections of unstable steady states (or specifically, unstable fixed points), while the black lines represent that for the asynchronous stable steady states (or stable fixed  points). We can also observe that with sufficiently low [DA], oscillations can disappear through a phase transition or bifurcation (specifically, a Hopf bifurcation) (Strogatz, 2001), such that the neurons in the network are tonically and asynchronously firing at stable rates. Clearly, we can see that higher  level can laterally shifts the onset of oscillation (bifurcation point) rightward (compare dotted and bold), which means a higher [DA] level is required to maintain the oscillations. Moreover, the range of oscillation amplitudes are also significantly more constrained. Globally, we can also map out the network behavior with respect to the [DA] and , i.e., a phase diagram. The phase diagram in Figure 5A2 clearly shows that oscillation behavior can occur   level. The inset in Figure 5A2 is a replicate of Figure 4B.  (Figures 5D1,D2), similar to Figure 5C. A similar minimal (maximal) oscillation frequency (amplitude) can be observed (Figure 5D2 inset, dotted). However, for high  ([5-HT] > 2 nM) and [DA] (>7 nM), oscillation behavior becomes more easily attainable (bold).
A similar analysis is done on Pyr3-type(D1+5-HT1A+5-HT2A) pyramidal cells, instead of Pyr1-type. Figure 6 summarizes the analysis of 2-population network of Pyr3 paired individually with the same four (Int1-4) types of inhibitory neurons. The phase diagrams (Figure 6, Figures 6C1,C2 (Pyr3 and Int3) seem to be a hybrid of Figure 6A (Pyr3 and Int1) and Figure 5C (Pyr1 and Int3), with one of the asynchronous regions ( Figure 6C2). This implies that when  is sufficiently high, oscillation can occur even with very low level of [DA] (<3.5 nM). The network with Pyr3 and Int4 ( Figure 6D) looks similar to that of Pyr1 and Int4 (Figure 5D), but the oscillation regime can now occur over a much larger range of [DA] and .
We have seen how [DA] and  can modulate the network consisting of different combinations of pyramidal cells and inhibitory neurons in Figures 4-6. Taken together, we can make several observations. Firstly, we can observe that pyramidal cells with (excitatory) 5-HT2A receptors can oscillate even with low [DA] (Figures 6A,C,D) as compared to pyramidal cells with no 5-HT2A receptors (Figures 5A,C,D). Secondly, inhibitory neurons with 5-HT2A receptors can enhance inhibition in the circuit, which can cause oscillation to cease (compare Figure 6D2 with Figure 5D2). Thirdly, higher [DA] will inhibit interneurons expressing D2 receptors due to the latter's inhibitory nature upon activation, and as a result, network oscillation will cease (Figures 5C2, D2, 6C2, D2).

HETEROGENEOUS NETWORK MODEL
After investigating the DA and 5-HT modulation on various possible two-population excitatory-inhibitory networks, we shall now study the neuromodulation and drug effects on the dynamics of a 8-population network fully connected with neurons expressing all the considered receptor types and their combinations.  (Figures 7A-C). This is especially pronounced for pyramidal cells which express D2 receptors (Pyr2, Pyr4, Int3, and Int4), exhibiting an inverted U-shaped modulation (Pyr2 in Figure 7B; Pyr4, Int3 and Int4 not shown due to the modulation on these neuron by DA is similar to that of Pyr2). The oscillation frequency decreases from low gamma to beta band with increasing [DA], before it slightly increases again with further increase in [DA] ( Figure 7D). The value of [DA] producing the minimal oscillation frequency (∼6 nM) coincides with that of the maximal neuronal firing rates.

5-HT and DA modulation
Next, we vary    Figure 8 shows that 5-HT modulates the network activity in an interesting manner. The network oscillates either at a low or high  level, while intermediate  level (within the 1.08-2.22 nM range) leads to asynchronous tonic stable activity (Figures 8A-C). The underlying reason for such a phenomenon is due to the different affinities of 5-HT1A and 5-HT2A receptors. This intermediate tonic stable state may not arise if In fact, we have observed such multiple oscillation regimes in the simpler two-population excitatory-inhibitory network model (Figures 5B2,C2, 6C2,D2). For pyramidal cells without 5-HT2A receptors (Pyr1 and Pyr2), its activity is almost fully suppressed by 5-HT1A receptor inhibition (Figure 8A; Pyr2 not shown). The slight increase in activity with oscillation is indirectly activated by other neuronal subgroups (e.g., excitation from oscillating Pyr3-type neurons; Figure 8B). The frequency of the oscillation increases with increasing  before the latter reaches a Hopf bifurcation point (1.08 nM), after which the oscillation ceases. When [5-HT] exceeds a larger critical value (2.22 nM), the frequency decreases toward a stable value of about 24 Hz ( Figure 8D). It should be noted that the activations of neurons at high levels of  are observed only for Pyr3-type (D1+5-HT1A+5-HT2A) pyramidal cells and Int2-type (D1+5-HT2A) inhibitory neurons, while the other types of neurons are inhibited (not shown). Finally, we vary [DA] and  simultaneously, and find that if [DA] is above a certain value (10.946 nM), the network always oscillates for any  level ( Figure 9A). Specifically, there exist a Λ-shaped green curve in Figure 9 below which the network cannot support oscillatory activity (black region), while above it oscillation occurs. Although the oscillation frequencies generally decreases with increasing [DA] levels, they stay around the same range (Figure 9). Moreover, the frequency of the oscillation non-monotonically depends on ; increasing before  exceeds the left branch of the Λ shape in the phase diagram, and then decreasing after  exceeds the right  branch of the Λ shape (Figure 9). When [DA] > 10.946 nM, there exists an optimal  value where the network can attain the maximum frequency oscillation. The peak of the frequency lies in the gamma range, which is often observed during attentional processing (Benchenane et al., 2011). Interestingly, when [DA] is smaller than 4.1 nM, there is only a narrow finite range of  values that supports oscillatory behavior.

5-HT and DA receptor selective agonist/antagonist
Having observed how the PFC network model can be modulated by [DA] and , we shall now investigate the influence of DA and 5-HT receptors selective agonist or antagonist. A selective agonist (antagonist) is a drug that can activate (block) a specific receptor without affecting other receptors. Our model can mimic the effect of receptor selective agonist or antagonist by decreasing or increasing the half maximal effective concentration of the receptors, namely, [DA] 1 for D1, [DA] 2 for D2, [5-HT] 1 for 5-HT1A, or [5-HT] 2 for 5-HT2A receptors (Lambert, 2004;Golan et al., 2007).
To investigate the effect of D1 selective agonist or antagonist, we fix [DA] 2 = 8 nM,    (Figures 10A,C) than that of D2 expressing neurons (Figure 10B), before they all reaches a stable asynchronous tonic state .
Next, we vary [DA] 2 while fixing [DA] 1 = 4 nM,  = 0.3 nM, [DA] = 4 nM, to mimic the influence of D2 selective agonist/antagonist on the oscillation. The results are shown in Figure 11. Unlike D1 modulation, no bifurcation happens when we vary [DA] 2 , and the behavior of the network does not change dramatically. The variation of [DA] 2 only slightly increase the oscillation frequency (Figure 11), but affects the oscillation amplitudes of only D2-expressing neurons (Figures 11B,D).
To study the effects of selective 5-HT 1A agonist/antagonist, we vary     2 from 0 to 5 nM, 5-HT2A receptors transit from the activated state to the blocked state. As a result, the amplitude of the oscillation decreases, but the frequency of the oscillation increases (Figure 13). It is to be noted that the different neuronal types are affected differently. For example, the firing rate of pyramidal cells expressing D1 and 5-HT1A receptors decreases to approximately 10 Hz (Figure 13A), but that of pyramidal cells expressing D2 and 5-HT1A (or 5-HT1A and 5-HT2A) receptors decreases to less than 1 Hz ( Figure 13B). The firing rate of interneurons expressing D1 and 5-HT1A (or 5-HT2A) decreases to approximately 20 Hz (Figure 13C), and that of interneurons expressing D2 and 5-HT1A (or 5-HT2A) receptors decreases to less than 5 Hz ( Figure 13D). As atypical antipsychotic drugs typically block 5-HT2A and D2 receptors (Maher et al., 2002(Maher et al., , 2011, we also investigated such drug effects by increasing the values of

SUMMARY OF RESULTS
In this work, we have shown, from single neuron to neuronal circuits, how DA and 5-HT, with their multiple receptors and combinations, can tonically modulate the PFC neural activity, resulting in a variety of complex behaviors.
Due to the different affinities and opposing effects of the 5-HT1A and 2A receptors, the neuronal firing activity of a PFC excitatory neuron coexpressing these two receptors can be inhibited before being enhanced as 5-HT concentration increases. When we extend our analysis to the two-population excitatoryinhibitory neuronal networks, we find that generally, pyramidal cells expressing D1 receptors can provide various interesting network behaviors. In particular, 5-HT and DA can modulate the amplitude and frequency of the network oscillations. Depending on the receptor types expressed by the neurons in the network, 5-HT and DA modulation can cause the oscillations to emerge or cease. Hence, this can result in a finite oscillatory regime, which can create optimal oscillation frequency and amplitude with   respect to certain [DA] or  level. Moreover, we find that certain combinations of receptors are conducive for the robustness of the oscillatory regime, and for the existence of multiple oscillatory regimes. The analysis of the two-population model provides us a leverage to understanding the more complex and realistic heterogenous network model. In the heterogeneous network model, the pyramidal cells and interneurons with all the considered combinations of D1, D2, 5-HT1A and 5-HT2A receptors are synaptically coupled. The model reveals that for the network to oscillate, it requires a sufficiently high level of [DA]. At intermediate levels of [DA], an interesting bimodal feature with respect to 5-HT concentration level appears -network oscillation can occur in two separate ranges of 5-HT concentration level. This bimodal feature is largely contributed by the different affinities and opposing effects of 5-HT1A and 5-HT2A receptors, as observed in the two-population model. Very low DA concentration level can suppress the oscillation regardless of the 5-HT level. Finally, we show that selective D1 receptor antagonists (agonists) tend to suppress (enhance) network oscillations, and shift from beta toward gamma band, while selective 5-HT1A antagonists (agonists) act in opposite ways. Selective D2 or 5-HT2A receptor antagonists can lead to decrease in oscillation amplitudes, but only 5-HT2A antagonists can increase the oscillation frequency.
Based on the analysis of the two-population and full network models, a general trend can be observed: the oscillation frequency will decrease if the change causes an overall increase in excitation within the network ([DA] 1 ↓, [DA] 2 ↑, [5-HT] 1 ↑, and [5-HT] 2 ↓), and vice versa.

RELATIONS TO NEUROPHARMACOLOGICAL DRUG EFFECTS
As mentioned earlier, abnormal beta and gamma band oscillations have been observed in various neurological and neuropsychiatric disorders (Spencer et al., 2003;Cho et al., 2006;Singer, 2006, 2010;Basar and Guntekin, 2008;Gonzalez-Burgos and Lewis, 2008;Gonzalez-Burgos et al., 2010). Schizophrenic patients (late responder to antipsychotic drugs) have been shown to have enhanced power in the beta2 (16.5-20 Hz) frequency band in the frontal cortex as compared to controls (Merlo et al., 1998;Venables et al., 2009). Using fluphenazine, an antagonist of both pre-and postsynaptic D2 receptors, beta2 in schizophrenic patients can be reduced (Kleinlogel et al., 1997). In our model, if we only simulate the effect of D2 postsynaptic receptor antagonist, the results is actually an enhancement of beta2 oscillation amplitude (Figure 11). However, antagonist of D2 pre-synaptic receptors can effectively decrease the overall DA concentration level, which can result in a shift in the oscillation frequency out of the beta2 range (Figure 7). Hence the model suggests that the antagonist effects on the D2 presynaptic receptors may be more dominant than the D2 postsynaptic receptors. In a rat model of Parkinson's disease, beta band oscillation in the frontal cortex is abnormally high compared to controls (Sharott et al., 2005). In that work, the authors showed that administration of apomorphine, a non-selective dopamine agonist which activates both D1-and D2-like receptors (but with higher preference for D2-like receptors), reduces the high beta band power and shifts the oscillation slightly toward higher frequency (from 28 Hz to about 35 Hz). In our model, D2 agonist generally reduces the oscillation amplitude of the beta band (Figure 11), while D1 agonist slightly decreases the oscillation frequency. The former is consistent with the experiment but not the latter. This discrepancy deserves to be further investigated.
Hallucinogens (psychedelics) are agonists of 5-HT2 receptors, enhancing PFC activity and metabolism in humans, and for treatment of psychiatric disorders (Nichols, 2004). It is shown that application of the 5-HT2 agonist (2,5-Dimethoxy-4-iodoamphetamine or DOI) as compared to 5-HT2 antagonist (ketanserin) can lower the power of the beta band in the EEG signal from the frontal cortex of anesthesized rats (Budzinska, 2009). Our model with high [5-HT] 2 value (mimicking 5-HT2A antagonists) also reduces the oscillation amplitude in the beta band (Figure 13), and thus is consistent with the experiments.

MODEL LIMITATIONS AND FUTURE WORK
Our modeling approach involves incorporating various experimental data to constrain the model parameters. This includes electrophysiological and pharmacological properties of PFC neurons and synapses, and how these are distinctly modulated by the different DA and 5-HT receptors. Thus, our approach lies more toward biologically constrained firing-rate models (Wong and Wang, 2006;Eckhoff et al., 2011) than abstract connectionist models (Fellous and Linster, 1988). This is a first step toward a systemic understanding of DA and 5-HT comodulation in the PFC. Admittedly, the model has its limitations.
As expected from large-scale biologically based modeling, many model parameters are involved here. We have tried to base as many parameters as possible from experimental data. Some of these parameters are directly obtained or inferred from experimental measurements, while others are based on indirect evidences or assumptions. The extensive investigations on the localization of DA and 5-HT receptors in PFC provided biological plausible proportions of subpopulations of neurons in the PFC network. However, we did not simulate all possible details in the model. For example, we did not include pyramidal cells which coexpress both D1 and D2 receptors. We also did not include PFC neurons which do not express D1, D2, 5-HT1A and 5-HT2A receptors. It remains unknown how these neurons will indirectly affect PFC network behavior upon 5-HT and DA co-modulation. Moreover, DA and 5-HT receptors generally have low and high affinity states, but the present model assume only receptors with high affinity states. In terms of selecting the parameters for the model, we have only chosen a single value for each parameter or variable (e.g., tonic basal  in the PFC) within a range of available values identified over various separate experiments. This problem is often encountered when integrating data from multiple sources during model development. Furthermore, oscillations in the cerebral cortex can differ among cortical layers, but the Frontiers in Integrative Neuroscience www.frontiersin.org August 2013 | Volume 7 | Article 54 | 15 current model does not deal with this issue. The present model also considers only tonic release state, which is more of a resting state and independent of any specific cognitive task. Despite these limitations and assumptions, it is sometimes advantageous to understand neuromodulation phenomena from simpler to more complex models, teasing apart the contributions of individual components of a system -a key advantage of computational modeling. As we can easily observe, even with such simplified models, the behaviors produced due to the DA-5-HT comodulation are already rather complex. Moreover, the specific [DA] 1 and [DA] 2 ,  1 and  2 values only reflect their comparative affinities with DA or 5-HT. Variation of these values while keeping their relative affinities will not dramatically change the network's qualitative behavior under DA and 5-HT co-modulation. That is, their absolute values are not as important as their relative values. A possible extension of our present work would be to explicitly specify the cortical layers, where the latter are known to be distinctively modulated by DA and 5-HT (Wang, 2010). Furthermore, for the model to generate slower oscillations such as theta and other lower frequency bands, and hence directly compare with other experimental data (Benchenane et al., 2010;Puig et al., 2010), the model may require additional slower dynamical features such as GABA Bmediated synaptic currents. These concerns will be addressed in future work.