A Neural Population Model Incorporating Dopaminergic Neurotransmission during Complex Voluntary Behaviors

Assessing brain activity during complex voluntary motor behaviors that require the recruitment of multiple neural sites is a field of active research. Our current knowledge is primarily based on human brain imaging studies that have clear limitations in terms of temporal and spatial resolution. We developed a physiologically informed non-linear multi-compartment stochastic neural model to simulate functional brain activity coupled with neurotransmitter release during complex voluntary behavior, such as speech production. Due to its state-dependent modulation of neural firing, dopaminergic neurotransmission plays a key role in the organization of functional brain circuits controlling speech and language and thus has been incorporated in our neural population model. A rigorous mathematical proof establishing existence and uniqueness of solutions to the proposed model as well as a computationally efficient strategy to numerically approximate these solutions are presented. Simulated brain activity during the resting state and sentence production was analyzed using functional network connectivity, and graph theoretical techniques were employed to highlight differences between the two conditions. We demonstrate that our model successfully reproduces characteristic changes seen in empirical data between the resting state and speech production, and dopaminergic neurotransmission evokes pronounced changes in modeled functional connectivity by acting on the underlying biological stochastic neural model. Specifically, model and data networks in both speech and rest conditions share task-specific network features: both the simulated and empirical functional connectivity networks show an increase in nodal influence and segregation in speech over the resting state. These commonalities confirm that dopamine is a key neuromodulator of the functional connectome of speech control. Based on reproducible characteristic aspects of empirical data, we suggest a number of extensions of the proposed methodology building upon the current model.


Introduction
Computational neuroscience takes a grounds-up approach to understand complex neural phenomena by investigating the underlying activity of neurons themselves. Starting from the basic neural model of Hodgkin and Huxley [1], which described the electric activity of a neuron in terms of the cell's constituent ionic currents, computational neuroscientists began to study temporal input-output relations of neural units. One of the most famous models created during this era was the ''point unit'' by McCulloch and Pitts [2]. Among the first modelers who incorporated not only temporal but also spatial aspects of neural processing was Rall, who used compartment models to show the strong impact of dendritic arborization on neural processing of synaptic inputs [3]. His work laid the ground for the first neural network modeling based on rather complex single neuron models. Next, the so-called Wilson-Cowan units [4] allowed for simulation of rather realistically macroscopic responses of entire brain regions on the scales corresponding to measurements obtained by non-invasive in vivo human imaging techniques. Indeed, neural simulations did fit well with data from various human imaging modalities, such as magnetoencephalography (MEG) [5], positron emission tomography (PET) [6], and functional magnetic resonance imaging (fMRI) [7,8]. Many neural states have been modeled, resulting in a rich literature on resting-state brain activity as well as behavior-specific activities, such as visual [6], memory [9], sensorimotor [10] and auditory [11] processing. However, while it is generally accepted that differences may be seen between the resting state and task conditions as well as between healthy and patient data and models, the fundamental question in modeling and data analysis methodology, the significance of differences between modeled conditions, remains unclear [12].
Given that functional activity may be affected by neurotransmitters (see [13] for a review), recent modeling efforts have been undertaken to integrate neuromodulators, such as dopamine, into task simulations. Dopaminergic neurotransmission has been implicated in cognition, learning, motor control, and, more generally, sensorimotor integration [14,15,16]. Chadderdon and Sporns proposed a large-scale computational model of prefrontal cortex to show the effects of dopamine release on the onset and performance of working memory tasks, which could be confirmed by behavioral, single-cell and neuroimaging data [17]. Determining how dopamine may regulate the functional connectivity observed during a behavioral task is a critical next step in addressing the ambiguity of task-specific functional connectivity.
To that end, we present a biologically-informed, large-scale model, which is based on neurobiological considerations, to simulate neuronal function and connectivity modulated by dopamine release in the human brain. The present paper applies this model to investigate speech production, one of the most complex human behaviors, which can be studied in a neuroimaging setting. Speech production is known to integrate several neural networks, ranging from auditory processing to motor control of articulatory movements [18,19]. Our recent study has demonstrated that dopaminergic modulation may play a role in left-hemispheric lateralization of functional brain activity and connectivity during speech production in the absence of lateralized structural networks [20]. Here, we propose an extension of a nonlinear model presented by Breakspear et al. [21] to allow for the simulation of brain activity due to dopaminergic modulation. We introduce the original model, which is a system of stochastic differential equations (SDEs), and rewrite it in terms of a multidimensional time-continuous stochastic process. Coupling between the brain regions with respect to regional neural firing rates is incorporated within the framework of Ito processes [ [22], Chap. 7]. A dopamine release model is developed and integrated into the model linking the basal ganglia and the laryngeal motor cortex based on previous studies [23]. We further present a mathematical proof establishing the existence and uniqueness of solutions to the extended model. Exploiting specific structural properties of the model, a computationally efficient scheme for numerical approximation of solutions is also presented. We show simulations of resting-state and dopamine-modulated BOLD signals and analyze the associated functional connectivity networks as related to corresponding real fMRI data obtained from healthy volunteers during the resting state and speech production. Finally, we discuss merits and limitations of the proposed model.

Ethics Statement
All participants provided written informed consent before participation in the study, which was approved by the Institutional Review Boards of the Icahn School of Medicine at Mount Sinai and National Institute of Neurological Disorders and Stroke, National Institutes of Health.

Modeling Objective
Our goal was to simulate a large-scale neural population using N coupled small scale local models, each replicating neural activation in a specific brain region i (i~1, . . . ,N), while incorporating neuromodulator release in a region-specific manner. Every regional subsystem consisted of interconnected excitatory and inhibitory neurons, which were assumed to be representatives of the local neural ensemble within a region. Thus, all quantities were understood as mean values across the considered region. The dynamics of regional state variables were governed by voltage-gated ion channels, functional synaptic couplings and neurotransmitter release. Thus, the temporal evolution of the entire population was determined solely by the interaction of its regional subsystems. In contrast to other approaches, the model discussed here was not based on coupled oscillator systems like the widely-used Kuramoto model (compare, e.g., [24] or [25]), but was based on neurobiological considerations. Below, we first detail the theoretical aspects of the model, including the Wilson-Cowan and dopamine dynamics, then describe the integration of the model with data.

The Breakspear Neural Model
Following Breakspear et al. [21], we denote the average membrane potential of neurons in region i by V i , which we assume to be governed by voltage-gated potassium (K), sodium (Na) and calcium (Ca) ion channels together with the passive conductance of leaky (L) ions. Thus, for j[ K,Na,Ca,L f glet m i j denote the fraction of open j-ion channels and let j be the ion population's maximum conductance for m i j~1 (i.e., when all j-ion channels are open). The basic model describing current flows across neural membranes in region i is a balance equation of the form (assuming unit neural capacitance) where V j denote the respective Nernst potentials and W j [C(R,R) are neural activation functions. To adequately reflect relaxation times of potassium channels, W K is characterized by an exponential decay with W 0 K being the value of W K at the initial time t~0, Q denoting a temperature scaling factor and t being the relaxation time. For brevity, we introduce the shorthand notation W i K :W K (m i K ). The other neural activation functions are defined as W CaW Na~i d(R) and W L :1.

Author Summary
Our knowledge of brain activity and network organization during complex motor behaviors in humans relies mainly on neuroimaging studies. However, the majority of available brain imaging methods are not feasible for quantifying the neural processes that occur on very short time-scales at the microscopic level. To address this shortcoming of functional MRI, we designed a mathematical model, which simulates brain activity using local ensembles of neurons and physiologically meaningful variables, such as cellular membrane potentials and ion channel relaxation times. We further incorporated dopaminergic function into our model as a neuromodulator of the dynamic organization of brain networks. We applied our model to examine brain networks controlling human speech and language production. We present a rigorous mathematical proof, which establishes the theoretical validity and solvability of the presented model, and we discuss the influence of dopaminergic transmission on simulated brain activity. We show that our model successfully reproduces characteristic changes seen in empirical data between the resting state and speech production. Our results indicate that the proposed mathematical model may be used as a platform for future studies to investigate the specific impact of certain pathologies within the dopaminergic pathways and their effect on global network dynamics.
Assuming that the ion channel specific opening-thresholds are normally distributed with mean T j and variance d j across the considered neural population, the fraction of open channels in region i may be computed as Note that the basic model (1) consists exclusively of uncoupled equations, i.e., the membrane potential of neurons in region i is entirely independent of neural firing in neighboring regions. Coupling is introduced by considering firing rates of excitatory and inhibitory neurons across the whole population. Thus let Z i be the mean membrane potential of inhibitory interneurons in region i, and define average excitatory and inhibitory firing rates as follows where V T and Z T denote the mean values and d V and d Z are the variances of membrane threshold potentials of excitatory and inhibitory neurons, respectively (assuming a normal distribution of thresholds across the neural population). Regional membrane potentials are altered by excitatory and inhibitory cell firing via synaptic feedback loops. Thus, functional synaptic factors a ee , a ie and a ei are introduced to scale excitatory-to-excitatory, inhibitoryto-excitatory and excitatory-to-inhibitory couplings, respectively. Furthermore, to reflect firing rate dependent glutamate neurotransmitter release (opening additional calcium ion channels) a supplemental scaling parameter r NMDA (the ratio of NMDA to AMPA receptors) is used. Non-specific input to excitatory and inhibitory neurons is modeled using random noise, which gives rise to a system of coupled stochastic differential equations (SDEs  (3) and (4) to be obviously adapted to V t (replace V i by the components V i t of V t in the respective definitions). In the following we establish a vectorial representation of the basic model equations given in [21]. Thus, for a vector u[R N let D(u)[R N|N denote a diagonal matrix with the components of u on its main diagonal. Further, we introduce the N-dimensional vectors and where a ni and a ne denote synaptic factors corresponding to nonspecific excitatory/inhibitory input, d is a noise scaling parameter which is the Ito version of the original multi-compartment neural dynamics model presented by Breakspear et al. [21].

Model Extension 1: Inter-Regional Connectivity
The original model (8) uses a scalar parameter C to parameterize excitatory coupling between regions. This means, in the framework considered here, that inter-regional connectivity strengths constant throughout the entire brain. To relax this restrictive assumption, we assign each pair of regions i,j f g coupling parameters c i,j and c j,i representing the connectivity strengths i?j and j?i, respectively. We collect the inter-regional coupling parameters in a N|N matrix C~c i,j È É N i,j~1 and incorporate it in the model (8) as follows. Instead of calculating excitatory-to-excitatory neural feedback by relying on a mean firing rate Q Q V , we scale neural firing using weight information from the coupling matrix C. Thus, we assume that firing of brain areas connected to region k impact the membrane potential in region k according to P N i~1 c k,i Q i V . Hence, (5) is modified to be We set c i,i~1 to reflect local excitatory input within a region. Note that we do not impose any restrictions on the directionality of regional couplings. Depending on the application considered, the above formulation allows the use of directed connections (i.e., a non-symmetric coupling matrix C) or undirected connections (C symmetric).

Model Extension 2: A Dopamine Release Model
A second extension to the original model was incorporated to simulate the effects of speech-induced dopamine release, as shown previously in real data [20]. We were especially interested in the effect of dopaminergic neurotransmission on the primary motor cortex [27] and its direct influence on the activity of the laryngeal motor cortex (LMC), which is a final common cortical pathway of speech control [28,29]. Keeping in mind the biologically-inspired channel model adopted in the present paper, elevated dopamine levels in the striatum (without a differential effect on either D1 or D2 type of dopamine receptors) were assumed to increase the probability that potassium, sodium and calcium channels of LMC neurons open, thus making these neurons more likely to fire. Hence, we simulated both D1-and D2-type modulatory effects on these channels [30].
We modeled the direct dopaminergic pathway from the substantia nigra, pars compacta (SNc) to the LMC [28]. Thus, we assumed that dopamine release was solely driven by neural activity in the SNc. Hence, let DA~DA t [R 2 Dt §0 È É be a two-dimensional Ito process with components DA ' t and DA r t denoting the dopamine concentration in the left (') and right (r) LMC respectively. We assume DA is governed by two simple mass balance equations To reflect the positive feedback of neural firing in the SNc on dopamine release we define where Q SN,h V denotes the neural firing rate in the left/right substantia nigra as defined in (4), and r : ½0,T? r min ,r max f g is a (timedefpendent) production rate. We assume that r attains a maximum value r max during speech production and is equal to a (positive) minimum value r min otherwise. The precise value of the uptake rate is taken to be a reasonable value from previous studies of extracellular dopamine levels [31]. Following [32], dopamine re-uptake was presumed to be governed by a Michaelis-Menten type kinetics equation where d max denotes the maximal uptake rate and k max is the Michaelis-Menten constant. Thus, a closed form representation of the considered dopamine model is As mentioned above, dopamine was assumed to affect the firing of the LMC by altering neural ion channel permeability. Thus, the effect of dopamine on potassium channels can be seen as a dependence of the gain in m LMC,h K on dopamine concentration. Hence we modify the equation governing the fraction of open potassium channels in the LMC as follows where b h t denotes a dopamine dependent gain. In the absence of dopamine we want a gain of unity, i.e., b h t~1 . Conversely, we also like to impose an upper bound on the gain. To achieve this, consider the expression where a[½0,1 is an antagonism parameter controlling the overall impact of dopamine on the gain b h t . Obviously, if DA h t~0 then b h t~blo , thus, by setting b lo~1 , a unity gain for a dopamine concentration of zero is established. Since 0ƒaƒ1 and assuming physiologically meaningful dopamine concentrations, i.e., DA h t ƒ1, b hi sets an upper bound for the gain. Finally, we modeled the impact of dopamine on calcium and sodium channels in the LMC using the gain b h t . We expressed the dopamine dependence of the permeability of those channels via varying the LMC's excitatory-to-excitatory functional synaptic coupling by introducing In the absence of dopamine we have b h t~blo and thus a LMC,h ee~a ee . Rising dopamine levels increase b h t and, in turn, a LMC,h ee until b h t reaches its previously established upper bound b hi , which gives a LMC,h ee~2 a ee . Thus, we have the estimate To establish a closed form representation of the full model, let for h[f',rg. Let further A 2 (t,V t ,Z t ,DA t ) be given by the right hand side of (6) and define Similarly, let B[R 2Nz2 be given by B~(B 1 ,B 2 ,B 3 ) T with B 1 and B 2 as defined in (7) and B 3~( 0,0) the full neural dynamics model can be written as where A is called drift (or deterministic force) and B is the diffusion (or random force) of the model. In the following section we discuss an efficient strategy to numerically approximate solutions to the model (20). A rigorous mathematical proof establishing existence and uniqueness of those solutions is presented later.

Numerical Approximation
We used time discrete approximation techniques to simulate sample paths of the SDE system (20). Extensive numerical experiments revealed pronounced non-linear dynamics of the model, which motivated the use of a higher order solution scheme. We encountered numerical instability of the widely used strong order 1.0 Milstein scheme [33]. Using a strong order 1.5 explicit Runge-Kutta (RK15) method, however, proved to be reliable. An explicit strong order 2.0 scheme yielded no notable improvements over the RK15 method but required a considerably higher computational effort. Thus, a RK15 scheme was specifically adapted to the model (20).
To establish a time discrete approximation of the solution to (20), we started by defining a discretization of the interval ½0,T. For K[N, let Dt~T=K be a step-size and define discrete time points t n~n Dt for n~1, . . . ,K. We introduce the Markov chain Y~Y n [R 2Nz2 Dn~0, . . . ,K È É to approximate the stochastic process X that satisfies (20). Thus we set Y 0~X IC and Dg n~gtnz1 {g tn . Note that (20) is a 2Nz2-dimensional nonautonomous SDE with constant additive scalar noise. This latter property is exploited to construct a highly efficient recursive solution scheme that has a considerably reduced computational cost compared to a general purpose SDE solver.
The following considerations are based on the family of solution schemes presented in [ [26], Sec. 11.2]. The vector form of an explicit order 1.5 strong scheme for a non-autonomous SDE with constant additive scalar noise is given by where and Dq n is a random variable representing the following double stochastic integral Rearranging terms, (21) can be simplified to Note that Dg n *N (0,Dt) and Dq n is also normally distributed satisfying as shown in [ [26], Chap. 10]. These properties play a key role in practice since they allow us to generate the pair of correlated random variables Dg n and Dq n in an efficient and straight forward manner: let r n and s n be independent N (0,1) distributed random variables, then Thus, an approximate solution of (20) was recursively computed following scheme (24) with auxiliary quantities (22) and (26).
Note that (24) requires three evaluations of the drift term A per step. In contrast, the Milstein method adapted to model (20) reduces to and thus requires only one function evaluation per step. However, unlike RK15, (27) reduces to an explicit Euler scheme in the absence of noise (zero diffusion). Thus numerical instability of the Milstein scheme for a model like (20) exhibiting pronounced nonlinear characteristics in the drift term was predictable. Note that it is possible to enforce convergence of (27) by substantially reducing the step-size Dt. However, this in turn dramatically increases the total number of time-steps making the overall computational performance of the Milstein method significantly worse than that of RK15 (24). Hence RK15 was the solver of choice for all simulations presented below.
In order to produce measureable changes in extracellular dopamine levels, which reflect rapid phasic dopamine release during a behavioral task or a pharmacological challenge, the dopaminergic axons must be stimulated at frequencies of 10-20 Hz or greater [34]. Because phasic dopamine release may reach high concentrations for brief periods due to concerted burst firing of dopamine neurons [34,35,36], we tested our model at a neural firing rate.20 Hz with different time-step sizes. We found that a small step-size of 0.1ms had the highest numerical robustness and showed the optimal temporal resolution of neural firing in order for dopamine release/re-uptake to set in gradually, without jumps.
The simulations shown below have been run on a Mid 2010 Mac Pro (262.66 GHz 6-Core Intel Xeon, 24GB DDR3 ECC RAM) under OS X 10.9.1. All codes have been written in Python [37] making extensive use of the packages NumPy, SciPy [38] and Matplotlib [39]. Computationally expensive sections of the code have been converted to C extensions using Cython [40].

Integration of Model and Data
Data acquisition. The raw model output was converted to a blood oxygen level-dependent (BOLD) signal and compared to functional brain activity data in healthy volunteers. We used fMRI data of 20 right-handed monolingual English speaking subjects with no history of neurological, psychiatric, voice, or respiratory problems (13 females, 7 males, age 53:2+10:1 years [mean+SD]) as reported earlier [20]. Right-handed volunteers were recruited in order to control for brain activity lateralization differences between right-and left-handed people. All scanning sessions were performed on a 3.0 Tesla GE scanner equipped with a quadrature birdcage radio frequency head coil (Milwaukee, WI). Data were acquired under two conditions: 1) a resting state, during which the subjects fixated on a cross, and 2) a task production, during which subjects were asked to produce meaningful, grammatically-correct, short sentences. Whole-brain resting-state (rs-fMRI) images were acquired using gradient-weighted echo planar imaging (EPI) (150 contiguous volumes with TR 2 s, TE 30 ms, FA 90 degrees, 33 sagittal slices, slice thickness 4 mm, matrix 64664 mm, FOV 240 mm, in-plane resolution 3.75 mm, duration 5 min). To assure (21) the resting condition, these images were acquired before the taskproduction fMRI within the same scanning session. Physiological recordings were carried out using a respiratory belt to measure respiration volume and a pulse oximeter to monitor heart rhythm and were sampled at 50 Hz with the recording onset triggered by the scanner's selection trigger pulse. For speech-production fMRI, whole brain images were acquired using gradient-weighted EPI pulse sequences (TE 30ms, TR 10.6 s (8.6 s task production, 2 s image acquisition), FA 90 degrees, FOV 2406240 mm, matrix 64664 mm, in-plane resolution 3.75 mm, 33 sagittal slices, slice thickness 4.0 mm) with BOLD contrast and a sparse-sampling event-related design. The subjects were instructed to produce short meaningful grammatically correct English sentences (e.g., ''We are always away'', ''Tom is in the army'') after listening to an auditory sample. The auditory stimuli were delivered within a 3.6 s-period and the subjects reproduced the sentences within 5 s, followed by a 2-s image acquisition. A total of 36 trials per task (i.e., sentences, resting fixation) were acquired over the five scanning sessions with the tasks pseudorandomized between sessions and participants. All fMRI data was pre-processed using AFNI software package [41]. For rs-fMRI, the anatomy-based correlation corrections (ANATICOR) model [42] was applied to remove hardware-related noise; respiratory and cardiac signals synchronized with the EPI data were used to remove physiological noise based on the retrospective image correction (RETRO-ICOR) model [43]. The resting-state residual time series were spatially smoothed by a 6-mm Gaussian kernel within the gray matter and normalized to the standard Talairach-Tournoux space. Task-production fMRI. For speech-production fMRI, the first two volumes were discarded, the EPI datasets were registered to the volume collected closest in time to the high-resolution anatomical scan using heptic polynomial interpolation, spatially smoothed with a 6-mm Gaussian filter, normalized to the percent signal change and the standard Talairach-Tournoux space. The task-related responses were analyzed using multiple linear regression with a single regressor for the task convolved with a canonical hemodynamic response function. Based on empirical studies [44,45], the whole brain was parcellated into 70 regions of interest (ROIs), including 64 cortical and 6 subcortical areas (Fig. 1A).
Coupling Matrix. The coupling matrix C was based on anatomical connectivity estimated from fiber tractography using diffusion weighted data from nine out of twenty healthy subjects described above. A single-shot spin-echo EPI sequence with TE 80 ms, TR 8.9 s, FOV 240 mm, matrix 120|118 mm, 68 contiguous axial slices, slice thickness 2 mm was used to acquire whole-brain diffusion-weighted images. A total of 60 noncollinear directions with a b-factor of 1,000 s/mm 2 were used to measure diffusion. One reference image was acquired with no diffusion gradients applied (b0 scan). Based on the same 70 ROIs, the DTI data were processed using the FATCAT Toolbox of AFNI software [46] following standard steps to construct an averaged 70|70 structural connectivity matrix. The matrix was normalized with respect to its largest row-sum and used as coupling matrix C for all simulations presented below.

Results
Below, we start by showing existence and uniqueness of solutions to the model (20) (Theorem 1). Once this fundamental result is established, we present simulations generated by the model and analyze it with respect to empirical fMRI data.

Existence and Uniqueness of Solutions to the Neural Dynamics Model
The following result guarantees unique solvability of the model (20).
Theorem 1. For m [ R and sw0 let X IC *N (m,s). Then the system (20) has a unique t-continuous solution X t .
Proof. If we show boundedness and Lipschitz continuity of A and B on R 2Nz2 , then existence and uniqueness of a solution to (20) follows from Theorem 5.2.1 in [22].
We start by proving that A and B are Lipschitz continuous. Obviously B as a constant trivially satisfies a Lipschitz condition. Since linear and trigonometric functions are differentiable (and thus Lipschitz continuous), we only have to show that the Michaelis-Menten kinetics equation (12) is Lipschitz continuous with respect to dopamine. Thus, for x,y[R a straight-forward calculation yields where we used the reverse triangle inequality and the fact that where c is a positive constant and D : D denotes some vector norm on R 2Nz2 : Since all norms on a finite dimensional linear space are equivalent, we prove (29) for the maximum norm D : D ? . We start by showing boundedness of all components of A 1 given by the right hand side of (9) with LMC equations (18). First, note that all firing rates (4) are bounded by Q max V~1 and Q max Z~1 respectively. Furthermore, the rates of open ion channels defined by (3) and (14) respectively are bounded by 1. Thus, the neural activation function for potassium channels (14) satisfies Weights connected to region i may be estimated by where C i denotes the i-th row of the matrix C. Thus, let T be a vector in R N and consider the i-th component A 1 i of A 1 as given by the right hand side of (9) with LMC components (18 (17), all components a i ee of a ee satisfy a i ee ƒ2a ee . Thus we obtain the following estimate for the first term of where we used (31) and the fact that m i Ca ƒ1. Note that all terms subsumed in the constant a 1 w0 are independent of t and j i . Similarly, we establish Using (30) we further obtain Finally, we establish and due to Q i Z ƒ1, for f~(f 1 , . . . ,f N ) T . Thus combining (33) -(37) yields where l [ R 2 . Analogously to (37) we compute and hence readily obtain Finally, by (12) we have DDA sink,h Dƒd max , and thus we get the following estimate for (13) for any where we used r max ƒ1 and Q i V ƒ1. Thus we obtain Combining estimates (38), (40) and (42) for J~(j,f,l) T hence yields This together with the definition (7) of B eventually gives which establishes (29) with c~ĉ czDBD ? and concludes the proof.
Having established existence and uniqueness of solutions to the model (20), we now present simulations corresponding to the resting state and dopamine modulation and compare them to empirical fMRI data.

Simulated Temporal Brain Dynamics
Using the coupling matrix described above, brain activity was simulated corresponding to the resting state and task-induced dopamine release. A list of all used parameters is provided in Table 1, which were taken from literature and scaled appropriately to reflect units used in this work or manually estimated based on previously published values [21,17,31]. Physiological variations across simulated brain regions were modeled by normally distributing inhibitory-to-excitatory and non-specific-to-excitatory synaptic coupling strengths using a fixed random number generator seed across simulations. This introduced the possibility of regionally desynchronized temporal dynamics in the model allowing simulated neural nodes to evolve non-identically over time in the absence of inter-regional coupling. Note that all simulations below were run with the same initial conditions and parameter values, i.e., starting values and parameters were identical for the resting state and dopamine-modulated speech-related simulations.
In both resting-state and task simulations, complex spatiotemporal patterns of activity emerged. V for fifty simulated speech cycles. Note that the propagation of firing rate changes acted as a neural feedback loop on the SNc itself in that repeated dopamine release caused different activity patterns than preceding cycles. In the task simulation, the LMC exhibited on average slightly higher firing rates than during rest (rest: 0:35+0:07, task: 0:36+0:06, compare also Fig. 1B) in agreement with the initial modeling assumption. To highlight that the proposed dopamine release model indeed shaped the dynamics of the entire neural population, the following section discusses changes in the correlative structure of simulated brain activity under dopamine modulation relative to the resting state.

A Simulated Functional Connectome
The raw model output was converted to BOLD signals as detailed above. Fig. 2 shows simulated and real BOLD signals for a selection of speech-related ROIs ( Fig. 2A,B). Simulated BOLD signals with and without dopamine modulation were compared to empirical resting-state and speech production fMRI data, respectively, in order to assess the global effects of dopamine modulation on the entire simulated neural population. To do so, we employed graph theory analysis to quantify variations in functional connectivity between the resting state and speech production. Thus, we first had to quantify statistical similarity between two time-series. We chose the normalized mutual information (NMI) [47] as statistical metric. Hence, for two random variables X and Y , let H(X ) and H(Y ) denote their respective Shannon entropies [48] and define where I(X ,Y ) denotes the raw mutual information between X and Y . Hence, unlike the original formulation of the mutual information I(X ,Y ), which is not bounded from above [ [49], Chap. 2], the NMI is normalized by the geometric mean of the entropies H(X ) and H(Y ). Thus, I I(X ,Y ) takes values between zero (two signals are independent) and one (two signals mutually depend on each other), permitting unambiguous comparison of values across data sets. Pairwise interactions in the simulated BOLD signals with and without dopamine modulation were quantified by computing NMI coefficients for each pair of ROI time-series. Analogously, NMI matrices were computed for the group-averaged resting-state and speech production BOLD data. This gave rise to four 70|70 NMImatrices (model rest, model speech, data rest, data speech) (Fig. 2C,D). Visual inspection of the matrices revealed larger variability in the model's correlative structure than in the corresponding empirical data. This might be partly explained by the fact that the empirical data were averaged across twenty subjects in an attempt to minimize subject-specific effects. Averaging a number of simulation runs would possibly decrease variability in the model; however, the aim of this study was to establish a qualitative assessment of the presented dopamine release model with respect to global effects seen in empirical data. In that respect, the proposed model, incorporating a single dopaminergic link between the SNc and laryngeal motor cortex, modulated neural activity of the whole brain to an extent that differences were observed between the structure of model's functional connectivity during dopamine release and the resting state. In addition, the model's prediction of empirical functional connectivity during speech production was in good alignment with the data.
In the following, we discuss simulated and empirical functional connectivity using the framework of graph theory. Interpreting functional connectivity matrices as graphs allowed us to not only reveal the functional topology and connectivity architecture of data and model but to also rigorously quantify the observed differences using well-established network metrics (see Supporting Information). By interpreting the 70 ROIs as nodes u i of a network with the associated NMI-coefficients representing the weights of the graph's edges, we constructed four weighted undirected graphs. Note that with the NMI being always nonnegative (contrary to the classical zero-lag Pearson correlation coefficient) a graph-theoretical analysis of NMI networks is straight-forward. Without the need to either extend classical metrics to negatively weighted graphs or consider negative and positive edges separately, most graph measures can be readily applied to NMI networks.

Graph Theoretical Analysis
The four weighted, undirected networks were analyzed following the concepts of functional integration, segregation, and influence [50]. As a measure of integration, we considered the local efficiency e i of a node u i , i~1, . . . ,N, quantifying a node's local communication performance in terms of inverse shortest path lengths within its neighborhood [51]. The degree of functional segregation was estimated using the weighted local clustering coefficient c i , which was calculated as the average geometric mean of edge weights in triangular motifs around u i [52]. Nodal influence was approximated based on nodal strength s i and nodal degree k i . A node's strength is the sum of attached edge weights, while its degree is defined as the number of connected edges. Clustering coefficient and efficiency were also compared to corresponding values of 100 conservatively-configured, null-model random networks. Normalized clustering coefficient c c i and efficiency e e i were computed by dividing c i and e i by the respective random network values. Statistical significance of differences in network metrics between the resting state and task production was determined using a paired two-sample permutation test at pv0:05 adjusted for family-wise errors (FWE) using the maximal statistic T max [53]. All graph metrics were calculated based on the full networks in their original density without applying any thresholding strategy. Since density-reduction techniques may severely deter network topology [54,55,56,57] and might thus dilute subtle differences between simulated and empirical functional connectivity patterns, the presented analysis is focused on the full unthresholded NMI networks as suggested by [58]. Graph metrics were computed using a Python port (pypi.python.org/pypi/bctpy) of the Brain Connectivity Toolbox for MATLAB [59].
Nodal influence. Figure 3 shows nodal strengths of the networks. We found a significant increase in strength when comparing resting state to task production in both data and model (both pv1e{4). While the simulated networks showed a higher average strength than the data in the resting state (model: 43:84+1:36; data: 38:88+0:52), the difference was less pronounced during task production (model: 57:12+0:89, data: 55:58+0:83). Examining the distribution of nodal strengths in the data, we observed a marked right-shift of the distribution during speech as compared to the resting state, clearly reflecting overall elevated strength in the speech production network. This right-shift was seen in the simulated networks too, although to a lesser extent. The data showed a narrower strength distribution than the model in the resting state, reflecting higher variability of NMI coefficients for the simulated BOLD signal without dopamine modulation (compare also the corresponding NMI matrices shown in Fig. 2C,D). Nodal degrees of the networks did not reveal any particular structure. All networks (model and data) were maximally connected, i.e., all nodes had maximum degree N{1, which means all pairwise NMI coefficients were non-zero. Note that, unlike Pearson's correlation coefficient (PCC), the mutual information does not only reflect linear correlation, but also dependencies in higher moments [60]. While a zero PCC only indicates that there is no linear relationship between the observed quantities, two time-series have to be approximately statistically independent for the NMI to be zero (compare, e.g., [61]). In other words, two signals have to show a stronger kind of independence to yield an NMI coefficient of zero. Given that fMRI-based functional networks are largely composed of highstrength nodes, are fully-connected, and may be indistinguishable from random networks if unthresholded [62], it was not surprising to see overall positive NMI coefficients for the data. It was also expected that simulated BOLD signals generated by a system of coupled but structurally identical equations show large NMI coefficients.
Network segregation. As mentioned above, the local clustering coefficient c i quantifies the average weight of connected neighbors of the node u i . The networks considered here had maximal connection density, i.e., each node was connected to all other nodes in the graph. In this case, c i is not influenced by the presence or absence of edges and is thus given by the geometric mean of N{1 edge weights adjacent to u i . Hence, the local clustering coefficient is solely dependent on the nodal strength. Thus, c i (Fig. 4A) exhibited qualitatively the same characteristics as s i (compare to Fig. 3B). In both data and model, we observed a significant increase in clustering during task production as compared to rest (pv1e{4) (data: rest: 0.56+0.01, speech: 0.81+0.01; model: rest: 0.63+0.01, speech: 0.83+0.01) Interestingly, compared to the data, the model showed on average higher values of c i in the resting-state simulation, while the dopamine-modulated run exhibited very similar clustering characteristics. To assess differences in network topologies in contrast to random graphs, we compared c i to the corresponding random network values and computed the normalized clustering coefficient c c i (Fig. 4C). We found c c i to be greater than one in the dopamine modulated simulation and the empirical speech production networks, while both data and model failed to show values larger than one during rest (data: rest: 0.81+0.01, speech: 1.16+0.01; model: rest: 0.91+0.02, speech: 1.19+0.01). This indicated an overall elevated segregation of simulated as well as empirical speech production networks in relation to random networks. Furthermore, for both simulated and empirical networks the difference in values of c c i between rest and task was found to be significant (pv0:0001). Interestingly, with and without dopamine modulation the model showed a very similar variability in both c i and c c i compared to the empirical networks. However, while the data exhibited similar peak frequencies during rest and speech, a decrease in the most prevalent values of both c i and c c i was found in the simulated networks (Fig. 4A,C). This was indicative of a slightly lower variability of c i and c c i in the dopamine modulated simulation.
Network integration. We considered functional integration of the NMI networks by evaluating values of local efficiency e i (Fig. 4B). For fully connected networks, e i , similar to the clustering coefficient, is completely determined by the nodal strength, since the shortest path between two nodes is always given by their connecting edge. Thus, e i showed similar characteristics as s i . We saw a statistically significant increase during task production as compared to the resting state (all pv1e{4) with the model showing on average higher values of e i (data: rest: 0.06+0.001, speech: 0.17+0.005; model: rest: 0.08+0.005, speech: 0.18+0.006). We normalized e i to analyze differences in network integration compared to a set of comparable random graphs. We found significant differences in local efficiency between both data and model during rest and speech (all pv1e{4). However, similar to the clustering coefficient, simulated as well as empirical networks showed on average a normalized efficiency smaller than one during the resting state, while data and model exhibited values larger than one during task production (data: rest: 0.51+0.01, speech: 1.50+0.04; model: rest: 0.73+0.05, speech: 1.62+0.05). Simulated and empirical networks showed a significant (pv0:001) increase in normalized efficiency when transitioning from rest to task, which was indicative of larger nodal integration in the networks.

Discussion
We presented an extension of a model of neural assemblies proposed by Breakspear et al. [21] to simulate dopamine release in the human brain during complex voluntary behaviors. In contrast to other large-scale neural modeling techniques based on coupled oscillator systems, our approach was grounded in neuroanatomy and physiology and thus allowed us to design a dopamine release model guided by biological considerations. We established unique solvability of the proposed model and demonstrated a computationally efficient strategy to numerically approximate its solutions.
In the context of the model, we assumed the difference between the resting state and speech production to be solely given by a modulation of dopamine levels in the LMC via a direct input from the SNc. Thus, the model was oblivious to task-related effects caused by any neurotransmitter other than dopamine. Importantly, we observed pronounced differences between the resting state and task production in simulations. This finding may be interpreted as an indication of the profound physiological impact of dopamine on brain dynamics.
It is remarkable that altered neural firing rates within the bilateral LMC only were sufficient for the entire simulated neural population to exhibit changes in its temporal dynamics. We attribute these observed task differences to dopamine driving neural dynamics via the coupling matrix C. Given real structural connectivity data as input, the strength of the model lies in its ability to reproduce observed properties of connectivity during a dopamine-modulated activity with a biophysical prescription for dopamine neurotransmission.
To quantify the impact of dopamine release on the entire neural mass, we used functional connectivity and graph theory analysis and interpreted the results as a stationary synopsis of the global effects of dopamine modulation. The graph theoretical analysis of the functional connectomes revealed a number of similarities between the model and data. Due to the slightly higher variability of NMI coefficients, especially during the rest, (Fig. 2) nodal strengths of the model without dopamine modulation showed more fluctuations than corresponding values for the data. Nonetheless, strength values of simulated and empirical networks were in good agreement with the model exhibiting slightly larger values. Clustering coefficients of both model and data also showed qualitatively similar attributes, when comparing the resting state to dopamine modulation. Thus, the model showed characteristics comparable to those of the data with respect to functional segregation and nodal influence.
Similarly, local efficiency showed good qualitative agreement between simulated and empirical connectomes, thus the model mimicked functional integration patterns seen in the experimental data. However, both model and data failed to show increased network segregation and integration compared to random graphs during rest but showed consistently larger values than null model networks during speech production (all e e i and c c i greater than one for i~1, . . . ,N). This may support earlier findings indicative of pronounced changes in network organization for speech control [63,64]. It should be emphasized that normalized efficiency is computed using the notion of shortest paths within a network. Due to the absence of zero-weighted edges in the considered networks, the shortest path between any two nodes in the graphs was given by their connecting edge, effectively side-stepping the notion of paths in a graph. A thresholding strategy to eliminate 'weak' edges (i.e., edges corresponding to small NMI coefficients) may have somewhat remedied this problem. It should be noted, however, that interpreting efficiency values obtained from functional networks that are based on statistical similarity between brain areas is not an immediately evident approach. Indeed, since NMI networks express not only direct but also all indirect couplings between regions, a path-based metric, like efficiency, may yield ambiguous results (see, e.g., [65,66]).
Nevertheless, decreasing connection densities in the networks would also yield non-trivial nodal degree distributions, opening another perspective on nodal influence within the networks. However, weight-based thresholding must be performed with considerable precautions so as to not deteriorate topological properties of a network. Since this work was mainly concerned with establishing a biologically-informed, large-scale model with optional dopamine neuromodulation, no thresholding strategy was applied to the constructed functional networks. A future study focused exclusively on the graph-theoretical analysis of functional networks should address this issue.

Limitations
A visual inspection of the simulated and empirical functional connectomes (Fig. 2) revealed that the model tended to slightly overestimate regional pairwise interaction during both resting state and dopamine modulation. This finding was not surprising considering the fact that the simulated BOLD signals were generated by 140 structurally identical equations that only differed in some parameter values. This is an apparent limitation of the presented approach. However, one of the advantages of the presented model is that the employed strategy enabled us to perform large-scale simulations of brain activity based on considerable neurobiological detail without becoming too complex to be practically unfeasible.
Coupling between regions in the model was achieved via scaling excitatory neural firing rates by entries of the coupling matrix C (compare eqn. (9)). Thus, all modeled axonal connections were excitatory, which is a simplification that ignores the effects of feedforward inhibition. In particular, firing of connected areas impacted a region's membrane potential through excitatory projections targeting local populations of NMDA and AMPA receptors. In other words, inter-regional coupling was not modeled as an explicit consequence of changes in neural voltages of neighboring areas. Instead, the influence of other regions on the local membrane potential was mediated by changes in neural firing rates. In this context it should be noted that the proposed model did not include an explicit representation of inter-regional axonal conduction delays. To some extent, however, the employed form of indirect coupling in the model may be interpreted as a lumped representation of conduction time delays.

Future Directions
Having tested the model for its efficacy in reproducing essential features of real data from healthy humans during speech production, the next step should be an examination of clinical relevance of the proposed neural population model. This may be achieved by incorporating 'lesions' into a simulated network of interest to investigate the extent of inter-regional influences coupled with dopaminergic transmission in a range of neurological and psychiatric disorders, such as Parkinson's disease, dystonia, schizophrenia, etc. Furthermore, the model's use is not limited to human applications [67] and may be applied equally well to animal models of disease and normal behavior, taking into account appropriate modifications for differences in animal and human dopaminergic innervation [68].
Since most parameters used here were taken from literature, some inferences about trajectories of isolated regions can be formulated based on the exhaustive analytical treatment of the original model [21]. In the original model, inter-regional coupling is introduced using a scalar parameter C that acts on the spatially Table 1. All parameters used in the model (including the neural and dopamine components) are provided according to their notation used in the paper, with their description, their value, and their basic units. averaged excitatory firing rates of all modeled nodes, i.e., C P N j~1 Q j V . The approach presented here uses not a scalar, C, but a matrix, C~fc i,j g N i,j~1 , to introduce coupling and thus expands the scaled mean field firing to be P N j~1 c i,j Q j V . This can be seen as a weighted average of firing rates. Thus, in the absence of dopamine and for a diagonal coupling matrix C the dynamic behavior of an isolated node can be reduced to the cases discussed by Breakspear et al. [21]. However, for a general non-diagonal coupling matrix, the dynamics become increasingly more complex. Moreover, as our simulation results indicated, dopamine modulation also had a pronounced impact on the overall behavior of the model. Thus, the extensions proposed here changed the dynamics of the original model in a non-trivial way. Thus, a rigorous dynamical analysis of the presented model would require a thorough study of the non-linear relation between DA h t and Q LMC,h V (h[f',rg) and an assessment of the influence of dopaminerelated parameter choices on the temporal evolution of the LMC nodes in terms of a full sensitivity analysis [69,70,71]. It was not within the scope of this work to present such an exhaustive analytical treatment of the model. Nonetheless, this poses an interesting direction for future studies.
Given the demonstrated differences in functional connectivity across the entire experimental time in simulations of resting versus speech conditions, the question arises as to what extent dopamine altered function on small versus long time scales within the tasks.
Our results indicate that dopamine may influence dynamics on long time scales. This may suggest that rapid temporal release of dopamine, as evidenced by the spontaneous dopamine release incorporated during each time-step in the model, may be involved in slow plastic responses. Thus, it is tempting to speculate that a future adaption of the proposed dopamine model might yield further insight into the learning and adaptation involved in voluntary behaviors, particularly given dopamine's involvement in learning and motivational behavior in other tasks.
Finally, models simpler than the one considered in this paper are capable of reproducing empirical functional connectivity. In fact, a recent study showed that a stationary model of resting-sate functional connectivity explains functional connectivity better than more complex models [72]. In modeling empirical functional connectivity as accurately as possible, the application of a complexity reduction technique [73] to the introduced highly non-linear model should be considered in order to derive a set of considerably simpler equations of statistical moments. On the other hand, it has been shown [73] that functional connectivity is essentially state-dependent and that local changes of activity in a set of cortical areas (due to external inputs, attention, neuromodulation, or learning) change the dynamical state of the brain network, thus modifying the correlations between the brain areas and introducing various levels of complexity. Along this line, while simpler models have a number of computational advantages (e.g., reduced computational load, easier estimation of parameters, simpler relationship between structure and function), their ability to simulate complex temporal activity patterns at various cognitive   scales (and in the context of simulated dopamine modulation) may be somewhat limited. This motivated the development of the proposed complex model to better understand empirical data and to make predictions about the different states of dopaminemodulated brain activity during voluntary behavior. Future work should be directed to a possible simplification of this model, while assuring its ability to accurately reproduce the complex biological patterns of voluntary behaviors.

Summary
We conclude that a regional model that includes dopamine release, reuptake, and modulation of ion channels significantly alters the behavior of an otherwise unmodulated, resting state neural population model. This work thus combines a small-scale basic cellular biology understanding of dopamine to alter macroscopic behavior of neuronal systems with nontrivial structural circuity, and presents meaningful global simulated fMRI network behavior. Region-specific analysis warrants the identification of specific effects of neuromodulation on task-based networks for speech and other dopamine-modulated voluntary behaviors.

Supporting Information
Text S1 Hemodynamic model and network metrics. Details regarding the hemodynamic model used to convert the raw model output to BOLD signals as well as a short description of the used network metrics and the construction of the random networks. (PDF)