An integrated modelling framework for neural circuits with multiple neuromodulators

Neuromodulators are endogenous neurochemicals that regulate biophysical and biochemical processes, which control brain function and behaviour, and are often the targets of neuropharmacological drugs. Neuromodulator effects are generally complex partly owing to the involvement of broad innervation, co-release of neuromodulators, complex intra- and extrasynaptic mechanism, existence of multiple receptor subtypes and high interconnectivity within the brain. In this work, we propose an efficient yet sufficiently realistic computational neural modelling framework to study some of these complex behaviours. Specifically, we propose a novel dynamical neural circuit model that integrates the effective neuromodulator-induced currents based on various experimental data (e.g. electrophysiology, neuropharmacology and voltammetry). The model can incorporate multiple interacting brain regions, including neuromodulator sources, simulate efficiently and easily extendable to large-scale brain models, e.g. for neuroimaging purposes. As an example, we model a network of mutually interacting neural populations in the lateral hypothalamus, dorsal raphe nucleus and locus coeruleus, which are major sources of neuromodulator orexin/hypocretin, serotonin and norepinephrine/noradrenaline, respectively, and which play significant roles in regulating many physiological functions. We demonstrate that such a model can provide predictions of systemic drug effects of the popular antidepressants (e.g. reuptake inhibitors), neuromodulator antagonists or their combinations. Finally, we developed user-friendly graphical user interface software for model simulation and visualization for both fundamental sciences and pharmacological studies.


Introduction
Neuronal activities, through the firing of action potentials and synaptic transmissions, can be modulated by endogenous neurochemicals called neuromodulators, acting through biophysical and biochemical processes [1,2]. These neuromodulators are released by a distinct population of neurons, and for example, by altering the affinity of the associated receptors that influences the release and reuptake mechanism of the monoaminergic systems [5,6]. As neuromodulators can also influence the biophysical properties of the neurons and synapses via multiple receptors with differential affinities, the complexity level in a neuronal circuit function can be substantial [7,8]. Experimental work often focuses on a specific brain region or system (e.g. certain receptor subtype) or employs a specific experimental methodology specific to the single level of biological organization (e.g. whole-cell recording at the neuronal level or voltammetric recording at specific brain region). Thus, it is difficult to reconcile their systemic implications.
Sufficiently realistic computational neural models can help us to integrate various data types from different studies, and can also generate testable predictions. However, modelling the detailed biophysical effects of neuromodulators can be complex and computationally costly [9 -11]. In particular, neuromodulation often involves intracellular signalling processes at the pre-synaptic and post-synaptic sites, and can subsequently affect neuronal firing activities [12,13]. Hence, such computational models that incorporate these biological processes can be time-consuming to develop, and with the multiple model parameters and equations, computationally intensive to evaluate while posing a considerable challenge in scalability.
In this work, we propose a novel neural circuit modelling framework to circumvent such difficulties. To develop a scalable model, we make use of neural population-averaged activity (mean-field like) or firing-rate-type models that describe how the neural population activities depend on the averaged effective neuromodulator-induced currents. The latter are determined by neuromodulators' concentration levels and the corresponding receptor affinities. Compared with other more abstract population-averaged firing-rate-type models [14][15][16][17], our model parameters describing the input-output functions and temporal dynamics are informed and constrained by data integrated from a variety of experiments, which include electrophysiology, neuropharmacology, radioimmunoassay, voltammetry and microdialysis.
We discuss such modelling approach in the context of developing and simulating a neural circuit model interconnected among the dorsal raphe nucleus (DRN), locus coeruleus (LC) and lateral hypothalamus areas (LHA), which are major sources of the important neuromodulators serotonin (5-HT), NE and orexin, respectively. The motivations for selecting these brain systems to model are that they are important in regulating physiological functions especially in arousal, are known to interact mutually with each other, and are the targets of several drugs [13,[18][19][20][21]. In particular, the neuropeptide Ox is known to play an important role in energy homeostasis, food intake and appetite regulation, neuroendocrine functions and sleep -wake regulation [22 -24]. The monoamine NE is suggested to be responsible for numerous functions, including stress response, attention, emotion, motivation, decision-making, learning, memory and regulation of sleep (e.g. REM) [25][26][27][28][29][30], whereas the monoamine 5-HT can affect several physiological functions that include eating behaviour, emotion, and sleep regulation [31][32][33]. Abnormal 5-HT or NE levels are implicated in mood disorders and anxiety [18,34].
The overlapping roles of these neuromodulators are not surprising, given their mutual interconnectivity, and any targeted neurons could themselves be sources of neuromodulators (figure 1). We have also purposefully selected the Ox system as a case study to demonstrate how we can model a neural system that may not be well characterized (when compared with NE).
A comprehensive, biologically faithful, yet efficient computational model at the neural circuit level would enable us to conveniently evaluate, account and predict, at a systems level, important measurable variables such as the concentration levels of the neuromodulators, the neural population firing rate activities, the effects of individual or combined drugs (e.g. reuptake inhibitors or antagonists) and their interdependencies.
The organization of the rest of the paper is as follows. We first describe the general modelling framework. Then, as an example, we demonstrate the steps in modelling three mutually interacting brain regions and discuss the simulation results including drug effects. Next, we describe our user-friendly software for simulating and visualizing the behaviour of such models. Finally, we summarize the results and discuss the implications of this study.

An integrated modelling framework
To develop a biologically compatible neural circuit model would require knowledge of electrophysiological properties of the composing neurons in a specific brain region, and the nature of the interactions among themselves and with other neuronal groups. This can include neurons which themselves release neuromodulators. Modelling the release-and-reuptake/decay dynamics of the extracellular concentrations of the neuromodulators would require information inferred from in vivo voltammetry or microdialysis studies at the targeted sites under neuronal stimulation. We would also need to know how the variation in neuromodulator concentration can in turn affect neural firing rate activities via neuromodulator-induced currents, hence requiring knowledge of firing rate -neuromodulator concentration or firing rate-current relationships (figure 2). These neuromodulatorinduced currents typically involve relatively slow metabotropic G-protein-coupled receptor (GPCR) types, e.g. G-protein-coupled inwardly rectifying potassium (GIRK) or transient receptor potential (TRP) type cation currents  Figure 1. LHA, DRN and LC interactions. Arrows: effective excitatory connections between any two areas; circles: inhibitory connections. Different colours denote different brain areas and their respective connection types and the targeted areas [27 -31,35]. [36,37], on targeted neurons, which alter the neuronal firing rate activities. As discussed, explicitly modelling such signalling pathway mechanisms can be complex and computationally intensive if large-scale neural circuits are involved.
To circumvent such challenges, we turn to phenomenological yet biologically faithful models to mimic the overall effects.
To begin the modelling process, we first model the neural activity for each brain region using neural populationaveraged activity [38]. Because the time constant of the typical neural population firing-rate dynamics is approximately 10-100 ms, it is much faster than the dynamics owing to neuromodulators, which is approximately seconds to minutes (table 1). Hence, we shall ignore the neural population dynamics and assume the system to be dominated by the slower neuromodulator-induced dynamics [57]. In general, different neuronal types can respond differently (in terms of firing rate activity) to the same current injection. In experiments, such relationship is demonstrated by the frequency-current ( f 2 I) relationship. In a similar vein, we can describe the ith neural population by the population firing rate-current ( f 2 I) curve or input-output function [58]: where f i is the population firing rate, F i is the input-output function and I total,i is the total averaged afferent current. Under typical physiological ranges, it suffices to use a threshold-linear function [58]: 0, and 0 otherwise. K i is the constant gain or slope of the input-output function, I 0,i is the threshold current for non-zero firing and I bias,i is the current coming from other brain areas. Thus, after a specific threshold value of the averaged afferent current, the neural population will be activated, and there is a linear relationship between the neural firing rate and the overall afferent current. We shall later show that this function fits the experimental data for Ox, 5-HT and NE neurons. In general, the afferent current I i can consist of several different types of currents mediated by the different modulators and their receptor subtypes. Each of these currents will be determined by the corresponding neuromodulator concentration levels and the receptor affinities (figure 2).
For example, suppose a neuromodulator y from region j induces a current I j!i on target region i, then we can describe the dynamics of the current by where [y] is the neuromodulator concentration level, and t j!i is the effective time constant owing to the applied neuromodulator which can be estimated from experiments. (If there are more than one receptor subtypes mediated by the same neuromodulator over the same brain regions, then we can specify the above variables further, e.g. by defining t j!i,R and I j!i,R for a receptor subtype R.) The value of t j!i can be deduced from the response dynamics of the induced current (or neural firing rates, if the induced current data are not available) upon infusion of specific neuromodulator at the targeted neurons. The input -output function G j!i can be described by the sigmoid-like function commonly used in pharmacology [59] where p j!i,LR and p j!i,UR determine the range of the neuromodulatory effect on the currents, and p j!i,LS and p j!i,S control the lateral shift and the slope of the neuromodulator response current function, respectively. The values of these parameters will be fitted to experimental data through firing rate-neuromodulator concentration or firing ratecurrent relationships. We used the standard nonlinear regression method, nlinfit from Matlab (The MathWorks Inc., Natick, MA, 2000). This approach particularly allows us to circumvent the complexity of actually simulating the intracellular signal transduction at the post-synaptic neurons. This post-synaptic current depends upon the extracellular neuromodulator release which is, in turn, dependent on the ( pre-synaptic) neural population firing rate of the source neurons. Thus, to close the loop in the model, we have to mathematically describe how the release-and-reuptake dynamics are affected by the neural firing rate of the neuromodulator source. We follow a mathematical form similar to that estimated from voltammetric measurements [40]: where ½y p,j!i is the per stimulus [ y] release (at the targeted area i from source j ). The rightmost term in equation (2.5) represents the reuptake rate, and is approximated from the Michaelis-Menten equation. Here, K m and V max are the Michaelis-Menten constants, with V max defined as the maximum uptake rate and K m is the substrate concentration where the uptake proceeds at half of the maximum rate. The value for both of these parameters can be obtained from experiments [40]. In voltammetry experiments, f j is typically an artificially applied high current stimulus frequency to stimulate the release of y. However, following our previous work [60], we can redefine it as the neural firing frequency of the neuromodulator source. Hence, the value of ½y p,j!i has to be adjusted from that in voltammetry experiments, and the exact value in the model can be obtained by constraining the overall basal activity levels of the system to be within the observed experimental ranges (see below). rsif.royalsocietypublishing.org J. R. Soc. Interface 14: 20160902 Table 1. Basal firing rate, neurotransmitter levels, dynamical time constants, and other model parameters for the LHA -DRN-LC circuits. Asterisk: [39], assuming V max , and per stimulus release at dorsal lateral geniculate (DLG) and LC will be same. Hash: [40] For the case of Ox, there is a lack of available experimental data, particularly its release and reuptake dynamics. Hence, we adopt the simplest mathematical form to describe the [Ox] dynamics [60], with only two parameters: where a is the rise factor and h is a constant decay rate, and both considered to be free parameters. The value of a is selected, so that the release of [Ox-A/B] at DRN or [Ox-A] at LC is close to the observed basal value (table 1). A summary of the general model construction process is summarized in figure 2. Such a modelling approach can allow multiple brain regions to be constructed, simulated and analysed. (See Methods for a simpler approach when only two brain regions are considered.) Overall, we have proposed an efficient and scalable approach by incorporating neuromodulator properties and dynamics into traditional firing-rate-type models. We shall next apply this approach to develop a neural circuit model involving multiple interacting neuromodulators.

An example with three interacting neuromodulators
We shall now demonstrate, as an example, the steps towards developing a neural circuit model of three interacting neuromodulator systems (lateral hypothalamus, DRN and LC) through three corresponding neuromodulators (Ox, 5-HT and NE), based on available experimental data and equations (2.1) -(2.6). These brain regions were chosen mainly because (i) they consist of different neuromodulator systems that can directly influence each other, (ii) they were targets of existing drugs, and (iii) we can demonstrate how one could model with incomplete information. It is important to bear in mind that different datasets from various experiments are used for model construction and validation. First, the frequency-current ( f 2 I) curves for neurons from the three brain regions are determined according to the available electrophysiological data ( f 2 I curves and typical baseline firing rate ranges). Using the threshold-linear function (equation (2.2)), the fitted parameter values for the LC's NE neuronal f 2 I curve are k LC ¼ 0:058 Hz pA À1 , I 0,LC ¼ 0:028 pA and I bias,LC ¼ 37:41 pA [46]; that for DRN's 5-HT neuron are k DRN ¼ 0:033 Hz pA À1 , I 0,DRN ¼ 0:13 pA and I bias, DRN ¼ 24:82 pA [10,44,45]; and that for LHA's Ox neuron are k LHA ¼ 0:2 Hz pA À1 , I 0,LHA ¼ 0 pA and I bias, LHA ¼ 11:5 pA [43].
Next, parameters of the induced currents and associated G functions are fitted to the experimental data through the concentration-response relationships for the change in firing rate induced by the neuromodulator. For example, the total afferent current induced by neuromodulators on DRN's neurons can be rewritten as a sum: The terms on the right are the currents owing to Ox-A/B from LHA and NE from LC (for simplicity, we ignore the current owing to autoreceptors and interneurons within each brain region, see [17,61]). The dynamics for each induced current are dependent on the concentration of the neuromodulator.
The Ox-induced current on DRN neurons (when NE-induced current on DRN is clamped) can be described by ). With these parameter values, we were able to obtain reasonable fits with respect to experiments for the various input-output functions.
We now proceed to model the effects of 5-HT on LC's NE neurons. However, there is a lack of experimental data on the direct effects of 5-HT on LC's NE neurons (experiments typically focused on how different 5-HT 2 receptor agonists affect the firing rate of LC neurons [63]). Thus, to estimate the 5-HT-dependent input -output function of firing frequency, we approximate the input -output function from other experimental data by restricting the basal activities to  0.11 f M and f LC 2.15 Hz (table 1). We considered the same sigmoidal shape (as defined in equation (2.4)) for all the neuromodulator-dependent firing-rate where q j!i,LR and q j!i,UR determine the range of the neuromodulatory effect on the firing rate, and q j!i,LS and q j!i,S control the lateral shift and the slope of the neuromodulator concentration-dependent firing-rate function H j!i , respectively. Considering these baseline values, the estimated parameter values are q LR ¼ 2:15 Hz, q UR ¼ À2:15 Hz, q LS ¼ 4:273 nM and q S ¼ 0:3 nM [42]. Similar to the 5-HT modulating effect on LC, the 5-HT response direct firing -frequency curve was not available for LHA neurons. Thus, we approximate the input -output function from other experimental data [37,47] and restricting the basal activities to [5-HT] 1.6 nM and f LHA 2:3 Hz. Then, the estimated parameter values are q LR ¼ 2:3 Hz, q UR ¼ À2:3 Hz, q LS ¼ 0:9 nM, q S ¼ 0:1 nM [36]. Then, we approximated the parameters for I DRN!LHA to be: p LR ¼ 0 pA, p UR ¼ 36 pA, p LS ¼ À1:55 nM and p S ¼ 0:4 nM.
With these values, we observe that there is no Ox neuronal firing even for 10 mM of [5-HT] (figure 3e). This is due to the strong inhibition caused by the induced inward GIRK current (approx. 32 pA), which eventually saturates (at approx. 35 pA) for higher  levels.
As the model has a hard threshold in the f-I curve, there is a sharp change within the 10-100 nM of . Thus, we could not obtain a perfect fit for the functions G DRN!LHA ð½5-HTÞ and hence the f LHA À ½5-HT curve. Similarly, the fitted parameters for the NE-induced GIRK currents on LHA, I LC!LHA , are p LR ¼ 0 pA, p UR ¼ 120 pA, p LS ¼ À5:39 nM, p S ¼ 0:4 nM and t LC!LHA ¼ 1 s [56]. We encounter a similar issue for higher [NE] level, i.e. no perfect fit for G LC!LHA ð½5-HTÞ (figure 3f ).
After determining the input -output functions and dynamics for all the currents, the final step is to integrate all three brain regions and their interactions. In general, the activities for the combined system will be different from the individual isolated systems. Thus, the baseline activities of the coupled system will be different from that observed from the individual systems. However, the remaining set of parameters, the neuromodulator release per stimulus frequencies, i.e. the [ y] p 's, can be adjusted to resolve this. We found that for values of ½5-HT p,LHA ,   After successfully constructing the LHA-DRN-LC circuit model, we shall demonstrate simulating neuropharmacological drug effects in the system. In particular, we shall focus on effects of Ox-1 receptor antagonist (SB-334867-A), SSRIs and/or SNRIs on the LHA-DRN-LC circuit.

Drug effect simulations
Pharmacologically, antagonists can be classified into two categories: competitive and irreversible antagonists [65]. Pre-treatment or application of competitive antagonist can shift the baseline dose -response curve horizontally. This shift towards the higher doses (of neurotransmitter) increases the effective dose (ED 50 ) value of the dose -response curve (where 50% of the maximal response of the dose is being observed). Conversely, application of an irreversible antagonist can cause shifts in the maximum range of the antagonist effect and does not affect the ED 50 value [66].
Ox-1 receptor antagonists have been suggested to encourage sleep, as well as treatment and prevention of many psychiatric disorders [67]. In particular, the Ox-1 receptor antagonist, SB-334867-A, acts as a competitive antagonist which rightward shifts the Ox-A response curve in 5-HT and NE neurons in DRN and LC [62,68]. Thus, we can easily incorporate the effect of SB-334867-A, by simply laterally shifting the function G LHA!DRN=LC ð½Ox À AÞ.
Selective serotonin/norepinephrine reuptake inhibitors (SSRIs/NRIs) are some of the commonly known pharmacological agents that are used for the treatment of various psychiatric disorders. The basic (acute) actions of these drugs are similar: primarily to increase the extracellular concentration level of their respective neuromodulator concentration by inhibiting the uptake process and reduce the synaptic clearance in the extracellular space [69]. There have been numerous studies conducted to understand the effects of uptake inhibitors on [5-HT] uptake [69,70]. In   [70] show that monoamine uptake inhibitors can affect the values of the Michaelis-Menten constants K m and V max in the limbic part of the brain. For example, 10 mM of fluoxetine, an SSRI (when applied to the substantia nigra area), increased the value of K m by about a factor of 5 but did not alter the value of V max significantly [69]. Thus, to incorporate the influence of SSRIs/NRIs in our model, we can mimic the different doses of SSRIs by varying the different values of K m . For simplicity, our model will ignore chronic or other long-term secondary actions such as receptor density changes. As we increase the value of K m,  (mimic SSRI), [5-HT] linearly increases in both the targeted areas LHA and LC (figure 5a,h solid blue). This increase in  level causes a significant decrease in f LHA (figure 5c, solid blue), which is consistent with experimental findings [36]. This in turn causes a decrease in [Ox-A/B] levels in the DRN and LC areas ( figure 5d,g solid blue). Interestingly, because of the network effect, there is a subsequent decrease in f DRN (figure 5f, solid blue), consistent with [71] even when we did not incorporate any inhibitory 5-HT autoreceptors [72,73]. However, 5-HT's effect on LC's NE neurons is minimal, consistent with [74], and therefore, f LC does not alter the [NE] levels in the DRN and LHA significantly (figure 5i, solid blue). These effects remained to be validated in the intact brain.
Next, we simulate the combined effects of SSRIs and NRIs, by increasing the value of K m, [NE] to three and five times its control value (400 nM) while varying K m,  as previously. The model shows that for higher values of K m,[NE] more  and [NE] are released in the targeted areas in the LHA and LC when compared with controls ( figure 5a,h dashed red and dotted-dashed pink). This suggests that other than K m,  , [NE] release in DRN also helps stimulate the release of  in these targeted areas, consistent with [75]. This rise in the  level significantly affects f LHA , ½Ox À A=B DRN , [Ox À A] LC , ½NE DRN and ½NE LHA while there is little impact on f LC ( figure 5b -e,g,i).
Finally, to assess the combined effect of SSRIs, NRIs and SB-334867-A, we set K m, [NE] to be five times the control value and mimic the influence of 10 mM SB-334867-A on DRN and LC (by changing the p LR value from 3.8 to 2 pA, p UR from 54 to 51 pA, shift factor p LS from 22.3 to 24.192 nM, and slope factor p S from 0.341 to 0.592 in LC). For DRN, p LS is changed from 22.08 to 22.97 nM and p S from 0.452 to 0.367. We find that this triple-drug combination can cause a further decrease in the f DRN and f LC , and a substantial reduction in [NE] DRN levels (figure 5e,f,i solid green), whereas the rest are not significantly affected by the addition of SB-334867-A (figure 5a-d,g,h solid green).

Software for model simulation and visualization
Using Matlab, we have designed and developed a software, called 'NModC' (neuromodulator circuit), with friendly graphical user interface (GUI) for simulation, analysis and visualization of the types of models described. The software is easy to use, and can easily be generalized to additional brain regions, other neural subpopulations and neuromodulator types. The user can visualize the activities of multiple brain regions dynamically and simultaneously. These brain regions are embedded in a rotatable three-dimensional glass brain using standard Montreal Neurological Institute (MNI) coordinates (figure 6a). The user can also further specify brain regions of interest to find the dynamical variables' time courses and mutual relationships (figure 6b) for more detailed analysis. The model parameters can be easily altered to visualize the variation in the steady-state values of the transients of neuromodulator concentration level and firing rates, and can also compare the firing rates of the two brain regions (figure 6c). Further details are described in the Methods section, and the software is available at https:// github.com/vyoussofzadeh/NModC.

Discussion
In this work, we have proposed a new computational modelling framework for incorporating essential biological features of neuromodulation in neural circuits. This provides a computational platform to link from low-level neurobiology to large-scale brain activities.
Our framework is based on the population-averaged firingrate type of model which has model parameters constrained by neurobiology. This is to be compared with other firingrate-type models without such constraints [14][15][16]. The model integrates pharmacological and electrophysiological data from separate experimental studies to constrain the input-output neuronal functions, and also the timescale and profile of the effective neuromodulator-induced currents.
Another key difference in our modelling approach is the consideration of the release-and-reuptake/decay dynamics of the extracellular neuromodulator concentration level that can be inferred from voltammetry. By doing this, we circumvented modelling the complex intracellular biochemical processes, but instead, directly modelled the concentrationdependent effect on neural firing rate activity based on pharmacological-electrophysiological data.
In particular, to allow our approach to be generalizable to multiple interacting brain regions, we introduced neuromodulator-induced currents to bridge the gaps between neuromodulator concentration levels and targeted neural firing rate activities-afferent influences from different brain regions can be accounted for by their summed currents. The timescale and response curve of the neuromodulator concentrationdependent currents were constrained by data from combined pharmacological and electrophysiological experiments.
We demonstrated our modelling approach with the example of a neural circuit that involves three mutually interacting brain regions (LHA, DRN and LC), which are also sources of three different neuromodulators: orexin, serotonin and NE, respectively. The transient and steady-state dynamics of the experimentally measurable variables (neural firing rates and neuromodulators' extracellular concentration levels) could be easily simulated. In particular, our model supported the coexistence of the observed (steady-state) baseline firing rates and neuromodulator levels found in separate experiments.
An important application of our model was the prediction of the effects of neuropharmacological drugs on neural circuits. We simulated the effects of SSRIs, NRIs and Ox-1 receptor antagonist on the LHA-DRN-LC model. We first showed that SSRIs could have a wide effect on the neural circuit, except the LC-NE system. Interestingly, SSRI could inhibit the DRN (decrease in f DRN ), the source of 5-HT, even though we did not implement its inhibitory autoreceptors. This effect was essentially owing to the direct effects on serotonin heteroreceptors on the LHA and LC, which in turn inhibited the DRN. Similar circuit-based effects could be explained for the addition of NRIs and/or Ox-1 receptor antagonists.
The constructed LHA-DRN -LC circuit model turned out to be dominated by a unidirectional influence between any pair of interacting brain areas (figure 5). Hence, these result in monotonic relationships (either increased or decreased) as the model parameters (e.g. K m 's) were varied. However, this need not generally be the case. For example, a more balanced (especially excitatory -inhibitory) coupled network could easily lead to emergent circuit oscillations or even non-monotonic effects [8,76]. In the latter case, it might then be possible to search for the optimal drug dosage. In fact, we had shown evidence of such nonlinear emergent behaviour when the model incorporated autoreceptors and non-principal (e.g. inhibitory GABAergic) interneurons (to mediate indirect connections) [17]. Moreover, the excitatory -inhibitory balance of the network can also be influenced by the co-release of the neurotransmitters (e.g. glutamate) [77,78]. In this case, our framework can still accommodate this by introducing additional dynamical equations to describe the effects of glutamate or GABA (for the same pre-synaptic firing rate).
Our work has also shown that administration of multiple drugs (serotonin/NE reuptake inhibitors and Ox-A antagonist) simultaneously can be simulated in neural circuits to search for the optimal mixture of drugs. However, the results remain to be validated as there is a lack of such work done in experiments. Hence, this will form model predictions that can help experimentalists in designing future studies. For example, multielectrode array in vivo recordings can be designed to study the wide-ranging effects of drugs on different brain areas. It would also be interesting to use the model to minimize the side effects of drugs, which is an important issue in neuropharmacology.
Our modelling framework is scalable to incorporate multiple brain regions and hence can be used to study large-scale brain effects. This includes studying the changes in REM/ non-REM stages or sleep-wake cycle [15,16], cortical dynamics [8] and cognitive-emotional processing [8,76]. This would require extending our current GUI software by including cortical brain structure and their connectivity with the neuromodulator sources. Thus, these models could potentially reveal insights into the relationships among various brain and behavioural disorders (depression, addiction, antidepressants and sleep disorders). Importantly, neuroimaging data, especially from positron emission tomography and functional magnetic resonance imaging, could potentially be incorporated into our modelling framework, bridging across multiple scales and modalities, similar in spirit to the popular dynamic causal modelling approach [79]. Interestingly, recent whole brain molecular imaging (functional magnetic resonance imaging (MRI)) of serotonin transporter to characterize 5-HT dynamics in humans before and after (e.g. SSRI) drug administration is now possible [80], opening up another possible application of our modelling framework.
In summary, we have proposed a promising new computational modelling framework that can integrate various experimental neurobiological data into a computationally efficient large-scale neural circuit model for simulating, testing and predicting the effects of multiple endogenous neuromodulators and neuropharmacological drugs.

A simpler modelling approach for modelling two brain regions
Note that if one is only interested in the mutual interactions of two brain regions, then one may ignore the induced current implementation step (figure 2, second column), and directly model the influence of [ y] on the firing rates f i [60], i.e.
where t j!i is the effective time constant owing to the injected neuromodulator y on the ith neural population. K i ([ y]) can follow a similar form as the G function in equation (2.4).

Numerical simulations
The neural circuit model simulation for the interaction of the three brain areas is computed by using the forward Euler numerical integration method which is applied to the set of the first-order differential equations. A time step of 1 ms was used. Smaller time steps were tested without affecting the results. These simulations can also be extended to other (e.g. second-order or fourth-order Runge -Kutta) numerical schemes. Simulations were run until stable steady states are obtained.

Graphical user interface software
To begin the software, the user presses the 'Start' button on the starting window of the GUI. This will result in the model outputs of the interacting brain regions (figure 6a). The outputs appear in the form of normalized neural (firing rate) at the specified brain regions, e.g. red colour represents relatively higher activity, whereas blue colour represents lower activity. The range of the colour map is based on the absolute range of [0 255] Hz. These colours of activity can change over time, reflecting their dynamics. The regions are embedded in locations based on the MNI coordinates in a glass brain. Three-dimensional rotation of the brain is also allowed in the software. Although the structure of the glass brain is currently based on normal human MRI data, it can be easily replaced by an animal (e.g. rodent) glass brain using the appropriate brain atlas. Once the model is converged after simulation, upon clicking on the 'Outputs' button in the starting window, the dynamical variables for the selected regions will appear in a new window (figure 6b). The variables are the (absolute) neural firing rates and neuromodulator concentrations of the selected brain regions. Both the individual variable's temporal dynamics and mutual relationships between the variables can be plotted. Upon clicking the 'Parameters' button in the starting window, a new window appears in which the model parameter values and the initial numerical values of the variables are shown. The model parameter values can be edited within this window. Once this is done, the user can resimulate the new model by pressing the 'Simulate' button within the same window. This generates the transients of baseline firing rates and concentration levels and shows the relationship among them (e.g. firing rates). For default values, all the transient activities eventually attain their stable steady states. These steady-state values of the system parameters vary as we change the model parameters. For example, to mimic the complex drug effects of SSRIs, K m is varied and corresponding changes in the steady-state values are analysed separately (see §2.3). If a mistake is made, then the user can always retrieve back the initial values of the parameters by clicking on the 'Default' button.