A Detailed Data-Driven Network Model of Prefrontal Cortex Reproduces Key Features of In Vivo Activity

The prefrontal cortex is centrally involved in a wide range of cognitive functions and their impairment in psychiatric disorders. Yet, the computational principles that govern the dynamics of prefrontal neural networks, and link their physiological, biochemical and anatomical properties to cognitive functions, are not well understood. Computational models can help to bridge the gap between these different levels of description, provided they are sufficiently constrained by experimental data and capable of predicting key properties of the intact cortex. Here, we present a detailed network model of the prefrontal cortex, based on a simple computationally efficient single neuron model (simpAdEx), with all parameters derived from in vitro electrophysiological and anatomical data. Without additional tuning, this model could be shown to quantitatively reproduce a wide range of measures from in vivo electrophysiological recordings, to a degree where simulated and experimentally observed activities were statistically indistinguishable. These measures include spike train statistics, membrane potential fluctuations, local field potentials, and the transmission of transient stimulus information across layers. We further demonstrate that model predictions are robust against moderate changes in key parameters, and that synaptic heterogeneity is a crucial ingredient to the quantitative reproduction of in vivo-like electrophysiological behavior. Thus, we have produced a physiologically highly valid, in a quantitative sense, yet computationally efficient PFC network model, which helped to identify key properties underlying spike time dynamics as observed in vivo, and can be harvested for in-depth investigation of the links between physiology and cognition.


Introduction
The prefrontal cortex (PFC) is a key structure in higher-level cognitive functions, including working memory, rule and concept representation and behavioral flexibility [1][2][3][4][5][6], and has been linked to impairments of these functions in psychiatric disorders like schizophrenia [7][8][9][10] or attention-deficit/hyperactivity disorder [11]. Our understanding of the computational and dynamic mechanisms underlying these cognitive functions, their neuromodulation, and their aberrations in psychiatric disorders, is still very limited, however.
Computational network models are a highly valuable tool for driving forward such an understanding, as data from many different levels of experimental analysis can be integrated into a coherent picture. With respect to psychiatric conditions, it is of particular importance that models incorporate sufficient biological detail and exhibit physiological validity in order to serve as explanatory tools. Psychiatric conditions like schizophrenia are characterized by a multitude of abnormalities in diverse cellular and synaptic properties, transmitter systems, and neuromodulatory input [7][8][9][10]. Moreover, pharmacological treatment options target the neurochemical and physiological level, yet they are supposed to change functionality at the behavioral and cognitive level. It is thus crucial to gain insight into the explanatory links between behavioral functions and the underlying neurobiological "hardware", a task that requires sufficient physiological detail in the model specification, in particular realistic assumptions about anatomical structure and cell type diversity.
Ultimately, the physiological validity of a computational model ought to be reflected in the degree to which it can reproduce and predict detailed aspects of the neural activity observed in vivo. That is, from a statistical perspective, one may define a good, physiologically valid model as one that accurately (i.e., quantitatively) captures distributions compiled from the electrophysiological activity (spiking, field potentials, membrane voltages) produced by networks in vivo, but not necessarily as one that captures every detail of membrane biophysics or receptor kinetics. In our perception, such requirements are currently not met even by sophisticated cortical network models which do include a lot of biophysical detail [12][13][14], as these are often only loosely compared to in vivo data or test only specific aspects of those.
In this work, we present a computational network model of the PFC which has high physiological validity and predictivity both at the single-neuron-(in vitro) and at the network-(in vivo) level, yet is still simple enough to be computationally tractable. Its anatomical structure, neural, and synaptic properties are completely derived from the experimental literature and our own experimental data. The activity of the network is compared with a range of statistics derived from in vivo data, including spike trains, local field potentials, and membrane potential fluctuations. The model turns out to reproduce these data quantitatively, and also exhibits robustness with respect to moderate changes in parameters.
in the real PFC (see Materials and Methods for details). The resulting model parameters (Table 1) follow broad distributions (Fig 1C), mostly of Gaussian shape, with the exception of Δ T and τ w which are best described by a Gamma distribution, and b which follows an exponential distribution (red curves in Fig 1C indicate distributions from which model parameters were drawn).
Anatomically, the network is divided into two laminar components, representing the superficial layers L2/3 and deep layer L5 (Fig 2A). Neurons are distributed over the five cell types in each layer based on estimates from the literature ( Table 2). The neurons are randomly connected with different connection probabilities p con for each pair of cell types according to the literature [18][19][20][21][22][23][24][25][26][27], including local clusters of higher connectivity [28,29]. The neurons are assumed to be organized in a single column and horizontal spatial distance is not taken into account. However, all neurons receive a constant background current (i.e., without fluctuations) that represents synaptic connections from outside the network, both within and outside the same column (see section "Admissable and realistic range of input currents" below). Since these currents were constant, all irregularity was produced intrinsically within the simulated network.
Neurons are connected by conductance-based synapses (AMPA, GABA A and NMDA) with kinetics estimated from electrophysiological data, short-term synaptic plasticity [30] that is matched to the types of the connected neurons [31,32], synaptic delays and random failure of synaptic transmission [33][34][35][36]. Distributions of synaptic weights (log-normal [37]) and delays (Gaussian) were extracted from the literature ( Table 3). The average connection strength (connectivity p con times synaptic peak conductance g max ) between pyramidal cells and interneurons in the different columns and layers is indicated by the width of the arrows in Fig 2A. Wherever possible, we used data from the rodent prefrontal cortex, or at least agranular cortices such as the motor cortex, which in rodents shows a similar layered anatomy as the PFC. Apart from the missing granular layer 4, specific features of the rodent PFC that are modeled here include an increased fraction of reciprocal compared to unidirectional connections [32], longer NMDA time constants than in other areas [38,39], and a uniquely prefrontal distribution of short-term synaptic plasticity properties for connections among pyramidal cells [32].

Reproduction of in vivo activity
To assess whether the network model can reproduce the dynamics of real prefrontal neurons in vivo, we compared measures computed from the model with those from electrophysiological  Arrows from or to one of the shaded blocks (rather than from or to a single neuron type) denote connection types that are identical for all excitatory (PC) or inhibitory (IN) neurons. Where all three types are drawn, they are randomly distributed over all synapses between these two neuron types according to the probabilities given in the figure. Right panel: Illustration of the postsynaptic potential in response to a series of presynaptic spikes for three types of shortterm synaptic plasticity for excitatory (E1 to E3) and inhibitory synapses (I1 to I3). data, as well as with a number of findings from the literature. Unless otherwise stated, we simulate a single column with 1000 neurons and apply a constant DC current of 250 pA to all pyramidal cells and 200 pA to all interneurons. These currents are the only parameters that are not directly obtained from experimental data. As discussed below, appropriate values for these currents were derived by inferring from lumped-population input simulations the amount of current produced by a network of realistic size, set up with the very same structure as the explicitly modeled network. Spike-train statistics. All experimentally recorded spike trains (kindly provided by Dr. Christopher Lapish, Indiana University Purdue University, Indianapolis, see [40] for details) were first segregated into statistically stationary segments to yield estimates of spike train statistics that reflect in vivo baseline activity, free from task-related responses (not modeled here) or Mean and standard deviation of the parameters of the synapses connecting the different pre-and postsynaptic neuron types (p con : connection probability, g max : peak conductance, τ D : transmission delay. other potential confounds [41]. For consistency, the same procedure was applied to the simulated spike trains, although, strictly, these were stationary by simulation setup. From all jointly stationary segments, the mean hISIi, coefficient of variation C V , and autocorrelation function of the inter-spike intervals (ISIs) were computed for each individual spike train, as well as the zero-lag cross-correlation CC(0) between pairs of neurons (Fig 3). The in vivo data show very low zero-lag cross-correlations between neuron pairs (2.4 Á 10 −4 ± 2.5, mean ± SD) and C V s near one (1.04 ± 0.33), consistent with the proposal of an "asynchronous-irregular" (AI) state of cortical dynamics (although the correlations theoretically proposed for the AI state are usually even at least one order of magnitude larger than obtained here [42]). The average single-cell ISIs follow a monotonically decreasing distribution with a mean comparable in size to the standard deviation (570 ± 610 ms), but with a heavy tail that is better described by a log-normal or beta-2 distribution [43] rather than an exponential distribution. The autocorrelation function shows a rapid decay with small negative flanks (halfwidth at half maximum: 10.1 ± 1.1 ms, minimum: 64.6 ± 69.9 ms, mean ± SD). Without further tuning of network parameters beyond their derivation from slice-physiological and anatomical data, all these in vivo statistics are well reproduced by the model (Fig 3). Two-sample Kolmogorov-Smirnov tests did not find notable differences between experimental and simulated distributions in any of the statistics (C V : p = 0.26, KS(29) = 0.28; mean ISI: p = 0.4, KS(29) = 0.23; CC: p = 0.4, KS(29) = 0.24), indicating that simulated distributions were not statistically distinguishable from the experimental ones. The asynchronous-irregular firing with low rates is also seen in the raster plot of spike times (Fig 3B).
Low fraction of spiking neurons and layer-dependent firing rates. Fig 3B reveals a relatively low fraction of spiking pyramidal cells in both layers-only 22% of the cells emitted more than 10 spikes during the 30s of simulated time, which will be used as the definition of "spiking neurons" throughout the paper, in line with [44][45][46]. Comparing the neural and synaptic parameters of those neurons which fire at a sufficiently high rate (> 0.33 Hz) and those which do not ( 0.33 Hz), we find that only the rheobase (and the cell parameters that contribute to it) differs between the two populations: Spiking neurons have rheobases at the lower end of the distribution (42.9 ± 2.1 pA), compared to 69.0 ± 1.6 pA for non-spiking neurons (mean ± SEM; p = 3.5 Á 10 −20 , t(997) = 9.4, two-sided t-test), some of them even firing spontaneously (called "generator neurons" [47]).
While neurons firing at very low rates may go undetected using extracellular single-unit recordings, recording techniques that are less biased toward spiking neurons, such as calcium imaging or in vivo patch-clamp, often reveal a large fraction of neurons that are mostly silent ("dark matter theory" of neuroscience, [44][45][46]). Consistent with these results, the fraction of neurons with more than 10 spikes rarely exceeded 40% in simulations with in vivo-like firing patterns (see section "Admissible and realistic range of input currents" below). This can be explained by the way the neurons are activated: While most neurons receive a background current above their rheobase, the high firing rates of the interneurons (Fig 3B) lead to an average membrane potential in the pyramidal cells below the firing threshold (mean difference: -17.3mV, range: -37.2 to -2.2mV for the example shown in Fig 3) that is occasionally kicked above threshold by random fluctuations. This means that the firing rate is mostly determined by the amplitude of the fluctuations of the membrane potential (see below for statistics). These results are qualitatively conserved across the range of input currents for which the overlap between experimental and simulated distributions is reasonably high.
Membrane potential and local field potential statistics. In addition to the spike data, we also compared the membrane potential statistics and LFP signals between simulation and experiments. For the simulated network, we observed a broad range of membrane potential fluctuations (after removing spike events; Fig 4A; 3.28 ± 0.72 mV, mean ± SD; range between 0.72 mV and 11.23 mV). We compared this distribution of standard deviations with those from in vivo patch-clamp recordings from 10 putative pyramidal cells during up-states in anesthetized adult rodent PFC (kindly provided by Dr. Thomas Hahn, Central Institute of Mental Health and BCCN Heidelberg-Mannheim). The simulated distribution is less than one SEM away from the average of the experimental distribution (pooled over all data sets) for most bins, and a Kolmogorov-Smirnov test (see Materials and Methods) does not show a significant difference (p = 0.45, KS(29) = 0.23). The range of membrane potential fluctuations in the model and in the recordings used here is also consistent with values found in the literature [48,49].
The local field potential (LFP) in the model was estimated as the sum of all synaptic currents (allowing excitatory and inhibitory currents to partially cancel). This is a reasonable approximation to the standard model of the LFP [50] under the assumption that all neurons are confined in a small volume of cortical space. We computed the power spectral density of this model-derived signal and of the LFP signals obtained from the in vivo recordings (Fig 4B). Up to a constant offset (that has been removed in the figure), the spectrum of the simulated LFP is less than one SEM away from the average estimated from the experimental recordings (from awake, behaving animals, also provided by Dr. Christopher Lapish [40]) at most of the frequencies. Both spectra follow a 1/f power law for frequencies below 60 Hz and change their scaling behavior for higher frequencies, consistent with LFP spectra described in the literature [51][52][53] (the fluctuations in the simulated curve are stochastic in nature, i.e. there is no systematic  deviation from the 1/f behavior across different simulations). For frequencies beyond 60 Hz, the experimental spectrum is well described by a 1/f 2 power law, while the simulated one rather follows a 1/f 3 relation. Both scaling behaviors have been reported in the literature (1/f 2 [52,54], 1/f 3 [51]), and the difference may result from the simplifications made in the computation of the simulated LFP, e.g. neglecting the spatial integration of currents in extracellular space or the contribution of active currents [14].
Transient information transfer and the role of neuronal heterogeneity. We next examined how neurons in L2/3 and L5 would respond to a simple stimulus simulated by a brief series of spikes at high rate (250 spikes within 5 ms) from a virtual (not explicitly simulated) "input population" connected to 10% of the pyramidal cells in L2/3 (cf. Table 3). The stimulus induces a number of spikes in L2/3, and with a short delay also in L5 ( Fig 5A). The delays (L2/3: 8.9 ± 1.1 ms; L5: 17.7 ± 1.2 ms, mean ± SD) are similar to values that have been reported in the literature (e.g. 3.4 ± 0.5 ms in L2/3 and 16.6 ± 1.2 ms in L5 [55]). Note that these delays are significantly longer than the fixed synaptic delays (below 2 ms, see Materials and Methods) and arise from the dynamics of the neurons and the kinetics of the synapses (c.f. [56]). For a sufficiently strong stimulus (e.g. 500 spikes within 5 ms), the neurons in L2/3 show a brief period (100-150 ms) of persistent activity ( Fig 5B).
The transmission of transient stimuli between layers crucially depends on the heterogeneity of the neuronal parameters. With a 80% reduction in the variance of all parameter distributions (but no change in the means), the stimulus only elicits a response in L2/3, but is not transmitted to the output layer L5 anymore ( Fig 5C). Indeed, L2/3 activity is almost independent of neuronal variability, whereas the number of spikes in L5 systematically decreases as the standard deviation of neuronal parameters is reduced (Fig 5D).
To further examine the transmission dynamics, we reproduced an in vitro experiment with suppressed inhibition [57] which showed that input in L2/3 resulted in an epileptiform spread of activation across the whole network under this condition, whereas the same input in L5 did not. We mimicked this setup by reducing the inhibitory synaptic weights in the network to 30% of their original values and inducing a strong stimulus (see above) in each of the two layers, while varying the peak conductance g max of the synaptic connection between the mimicked Poisson input population and the network. For moderate connection strengths (g max = 2), only a fraction of the network responds, and the number of spikes elicited by the network is much larger if the stimulus is injected in L2/3 (404 ± 116, mean ± SD) compared to a stimulus in L5 (118 ± 33). Higher connection strengths (g max = 20) reliably drive the network into an "epileptic state" (transient high-rate response from all neurons in the network) for a stimulus in L2/3. In contrast, this state was never reached for an input in L5, consistent with the experimental results in [57].

Conditions for in vivo-like dynamics
In the previous section we showed that the model can reproduce a wide range of characteristics of neural activity in vivo. Here, we assess how the reproduction quality of in vivo-like behavior depends on those parameters of the model which were only loosely constrained by experimental data. We restrict this analysis to the spike series statistics hISIi, C V and CC(0).
Admissible and realistic range of input currents. The background currents I ¼ ½I L23 ex ; I L5 ex ; I L23 inh ; I L5 inh have so far been treated as free parameters, as such estimates are difficult to obtain or at least have not been reported experimentally. We address this in two ways: First, we systematically vary these four currents and assess the similarity between experimental and simulated spike time distributions using Kolmogorov-Smirnov statistics as before. Second, we estimate the required background currents from the simulation itself, using the assumption that the simulated network is embedded in a larger, but structurally identical network from which these currents originate. Fig 6A shows the Kolmorgorv-Smirnov test statistic D KS as a function of I ex and I inh , where see below for a discussion of laminar differences in the input currents). The figure reveals that the overlap between experimental and simulated distributions is acceptable (p > 0.05 for the two-sample Kolmogorov-Smirnov test, i.e. failure to reject the null hypothesis H 0 of equal distributions for the two samples, see Materials and Methods) for a wide region of I ex and I inh values (delimited by the black isocline in Fig 6A, associated with D KS values below 0.4). More specifically, simulated C V and mean ISI distributions become indistinguishable from their experimental counterparts as I ex increases, while the overlap with the ISI distribution decreases again for very high I ex values (Fig 6A, left inset). Both C V and mean ISI deviate from the experimental distributions as I inh increases. CC, on the other hand, matches well with the experiments for high I inh values (Fig 6A, lower inset). As mentioned above, the fraction of firing neurons is quite low in most networks showing in vivolike firing patterns, typically between 20 and 30%, as shown in Fig 6B (blackly delimited region gives the empirically acceptable parameter regime copied from Fig 6A). The ratio of inputs into the two layers, I L23 ex =I L5 ex and I L23 inh =I L5 inh , does not have a strong influence on these results within the tested range (mean D KS ± SEM: 0.31 ± 0.05, 0.32 ± 0.05 and 0.35 ± 0.06 for ratios of 1, 2 and 4, respectively), but does of course affect the relative firing rate between the two layers. In vivo experiments found that firing rates are considerably higher in L5 compared to L2/3 pyramidal cells (3-20 times [27]). This condition is fulfilled in our model as long as L2/3 receives less or the same input as L5 (Fig 6C).
To estimate which range of I values could be realistically assumed, we tested whether a substantially larger network than the 1000-neuron-network simulated here would produce mean synaptic currents that are large enough to self-sustain in vivo-like activity (i.e. within the blackly circumscribed regions in Fig 6A). In this case, the activity in the large network and the small network (the latter driven by the currents from the larger one) would be indistinguishable, and the in vivo-like activity would be supported by the larger network. We increased the size of the network either by changing the density of neurons or by adding input from nearby columns (see "Estimation of background currents" in Materials and Methods). . The shaded areas show the ranges for I ex (blue) and I inh (red) within which these currents would produce in vivo-like activity (D KS < 0.4). Note that it is sufficient that one of the two layers receives a current above the lower bound, as it will push the other layer into the right regime by cross-layer synaptic connections. The upper bound, on the other hand, may not be exceeded by either of the two layers, as this would push the other layer beyond its upper bound as well. It is apparent that these conditions are fulfilled already for (spatially) relatively small networks (* 5 columns), and currents saturate as network size grows further (Fig 7A). By increasing the neuron density, on the other hand, the input currents increase monotonically over a wide range (Fig 7B, averaged over all column numbers ! 5). Mean synaptic currents sufficient to drive the network into the experimentally observed regime arise for densities between 19,000 and 44,000 neurons per mm 3 . This range overlaps with densities found in anatomical studies (30,000 to 90,000 neurons per mm 3 [58][59][60]; horizontal dotted lines).
Variation of synaptic parameters. We attempted to estimate all synaptic parameters from data reported in the literature. Given that these come with some uncertainty and variation, however, we explored how sensitive the network behavior is with respect to changes in mean synaptic peak conductances and their distribution, synaptic time constants, and the GABA A reversal potential. All these parameter variations were performed for a range of different background currents and averaged results are reported.
The GABA A reversal potential E GABA rev was initially set to -70mV, which is well within the range of the values reported in the literature [19,24,55,61]. Within the physiologically reasonable range from -90 to -60mV [62], the divergence between simulated and experimental  (Fig 8A). At the same time, the standard deviation of the membrane potential decreases. The time constants of the synaptic kinetics also turned out to be important for the agreement with in vivo data: While small changes are acceptable, both very fast and very slow GABA A kinetics strongly diminish the agreement with the experimental data (D KS = 0.99 for τ on = 0.6 ms and τ off = 8 ms and D KS = 0.85 for τ on = 12 ms and τ off = 160 ms). The NMDA time constants have less effect, unless they are very strongly increased (D KS = 1.0 for τ on = 17.2 ms and τ off = 300 ms, compared to values of τ on 7 ms and τ off 100 ms reported in the literature [38,39]). The effects of the   Fig 8B. While small to moderate changes ( ± 50%) have no significant effect, a strong decrease in the inhibitory synaptic efficiencies leads to a significant mismatch with the in vivo statistics.
Apart from the mean, we also analyzed how the distribution of the synaptic peak conductances affected in vivo-like behavior by either reducing the variability or drawing them from a normal rather than a log-normal distribution (conserving mean and standard deviation). Reducing the variability of the synaptic weights increased the mismatch between empirical and simulated distributions (Fig 8C, gray line). Surprisingly, just changing the form of the underlying distribution from log-normal to normal, without changing its mean or standard deviation, had a similarly strong effect as a pronounced reduction in standard deviation (Fig 8C, black  line), so both the variability as well as the functional form of the synaptic conductance distribution are crucial for reproducing spiking dynamics as observed in vivo. Without the long tail of the log-normal distribution, the network activity becomes much more synchronized (CC(0): 0.017 ± 0.006, mean ± SD) and exhibits strong bursts (C V : 1.81 ± 0.09, mean ± SD), while mean firing rates are not much affected (hISIi: 420 ± 558 ms, mean ± SD).

Discussion
We presented a model of the prefrontal cortex which is entirely defined by electrophysiological and anatomical data, and is capable of reproducing a wide range of in vivo statistics, including properties of single spike trains and pairwise correlations, the power spectrum of the local field potential and the variability of the membrane potential. Importantly, this reproduction did not require specific tuning of model parameters towards in vivo behavior. In fact, variation of the synaptic parameters shows that the ability of the network to show in vivo-like behavior is robust against considerable changes. In particular, it is possible to increase or decrease the synaptic weights by up to 50% of their original value without significant changes in the prediction performance ( Fig 8B). This keeps the model flexible to synaptic plasticity, i.e. the weights can be modified by task-related learning rules without changing the global activity of the network. The only variable that was not tightly constrained by in vitro data, the external input currents, showed a wide range of admissible values, within a range that would be produced by a network supposedly large enough to self-sustain activity. Despite the high biological validity, the model remains simple enough to allow for efficient simulation, due to the combination of a simple, but versatile neuron model and a complex network structure. To our knowledge, no other biophysical model currently exists that is as tightly constrained by the specific neuronal and synaptic properties of the prefrontal cortex and systematically compared to in vivo data as the one presented here.

Relation to other network models
The current model has a strong focus on its tight connection to data. Many existing network models of the neocortex are based on neurobiological findings as well [13,14,[63][64][65], but the present model differs from them in two respects: The strict way in which the in vitro data is used to fix or systematically infer every detail of the model, and, more importantly, the observed standard deviations (cf. Fig 4A). quantitative test of the model's validity on a wide range of in vivo findings. Recently, a few studies have also moved in this direction. Fisher et al. [66] proposed a model for the short-term memory circuit in the oculo-motor system of the adult goldfish. They fitted the model simultaneously to a range of anatomical, physiological and behavioral data. This approach gives a coherent picture of this particularly well-defined non-cortical system. Furthermore, Potjans and Diesmann [27] proposed a model of a sensory cortex network where the connection probabilities are thoroughly derived from in vitro studies. While the neuron parameters are generic and homogeneous, their focus is on the precise laminar and horizontal organization of the synaptic connections. They compare the results of their simulations with the baseline firing rates and flow of transient information through the different layers in vivo. These comparisons to experimental data are qualitative in character, as it is the case for most existing large-scale models of cortical networks [12][13][14]67]. However, a few recent studies also made statistical comparisons on partial aspects of physiological data [68,69]. It would be interesting to assess these models on a wider range of in vivo data as we proposed here, to see which degree of biological detail is sufficient to predict their key properties.

(B) D KS values as a function of percent change in overall synaptic peak conductances between pyramidal cells (E) and interneurons (I). The dotted line denotes the critical D
An important simplification made in the present model is the reduction to two laminar components, leaving out layer 4 and 6 as well as the long-range fiber bundles and interneurons in layer 1. While layer 4 is missing in rodent PFC, layer 6 is only weakly connected to the other layers in our reference connectivity maps, which are based on the motor cortex [57,70,71]. Thus, its inclusion in the network should not have a major impact on the results shown here. This is probably different in sensory networks, where layer 6 strongly interacts with both pyramidal cells and interneurons in layer 4 [72].

The relevance of synaptic and cellular heterogeneity
The model exhibits a low fraction of spiking neurons, consistent with results from recording methods such as calcium imaging, which are not biased towards high firing rates ("dark matter theory" of neuroscience [44][45][46]). As described above, this may partly result from the variancedriven firing of the neurons: The membrane potential is on average well below the spiking threshold, but large fluctuations still lead to occasional spiking. The size of the fluctuations and the low-rate, Poisson-like firing (C V % 1) of the neurons is consistent with the high-conductance state theory [48] and balanced-state theory [42,73]. We note that the irregular and highly asynchronous firing of the neurons [74] observed here is a generic property of the network that simply emerged from its parametrization through in vitro and anatomical findings.
There are two main determinants of the high-fluctuation regime of the model: First, variability in the membrane potential requires variability in the synaptic parameters and in particular, the fat tail of the log-normal distribution of the synaptic weights. Second, the range between the firing threshold V up and the GABA A reversal potential E GABA rev must be sufficiently large, because below E GABA rev , all synaptic currents depolarize the cell, so the dynamical range for a balanced, variance-driven state is constrained between these two values.
Using the multivariate distributions of neuron parameters obtained from our in vitro recordings, we also observed that decreased cellular heterogeneity has a profound effect on the processing of transient stimuli. It prevents the transmission of stimulus-induced activity from L2/3 to L5. This phenomenon can be understood if one considers the rheobase distribution: Reduced heterogeneity removes those neurons that originally had a very low or even negative rheobase. These are the ones which are highly susceptible to even small inputs and form a small but significant fraction of L5 neurons that were activated by the transient synaptic input from the L2/3 cells. Given that L5 provides the majority of output to other brain areas, impaired transfer of stimuli to this layer may lead to major impairments in information processing.
Thus, apparently quite subtle changes in the distributional properties of synaptic and cellular parameters (not affecting their means) may lead to major changes in network dynamics and functional connectivity among columns or areas, effects that have been proposed to underlie major psychiatric conditions like schizophrenia [8,9].
Self-sustained activity of the PFC By varying the total input from a virtual population designed according to the same principles as the actually simulated network, we provided evidence that a larger network than the one actually simulated with anatomically realistic neuron densities should be capable of self-sustaining in vivo-like spiking modes. Although we did not demonstrate self-consistency in a strict sense, we have shown that the background currents into the smaller, simulated network needed to yield in vivo-like behavior are consistent with the range produced by a much larger network of anatomically reasonable size. For currents within the blackly delimited region of Fig 6A, the spike train distributions are statistically indistinguishable from the in vivo statistics, and the background currents that would result from scaling up the simulated network to anatomically realistic size lie exactly within this regime. This analysis implies that in vivo-like activity can be self-sustained in a larger network with the same anatomical layout as explicitly simulated here, as it has been observed for instance in deafferented cortical slabs [75], while e.g. the thalamus or other sub-cortical structures may provide transient, stimulus-related input or modulate the overall activity of the network [76].
Interestingly, the currents produced by this procedure are much higher in L5 compared to L2/3 (Fig 7), as required for the much higher firing rates observed in vivo [27]. At first glance, this seems counterintuitive, as L2/3 neurons receive input from neighboring columns, while L5 neurons do not (see "Estimation of background currents" in Materials and Methods). However, L5 also receives strong inputs from L2/3, while the inverse projections are much weaker (Fig  2A). Thus, once L2/3 neurons receive enough input from other columns to spike, they drive L5 much stronger than themselves.
In terms of space, input from just a few columns is sufficient to drive the network, as connectivity rapidly decays over the cortical extent. Nevertheless, a single column is not sufficient for driving the network because of the higher fraction of excitatory synapses in long-range connections and the more local connectivity of interneurons. This is consistent with recent experimental studies [60,77] and earlier results from deafferented cortical slabs [75] (but see [78]).

Possible applications
In this study, we have focused on the resting state of the network. However, it may also be used as a foundation for more functional investigations of cognition. For instance, the clusters of increased synaptic connectivity may serve as building blocks for cell assemblies [29] which can be used to represent behavioral rules [3,40] or transient stimuli that need to be kept in working memory [64]. Moreover, the fast and fully automatized framework for fitting the neuron model to in vitro data [15] opens a convenient way to test the network effects of genetic or pharmacological manipulations: Recordings from neurons that underwent such a manipulation can be used by the very same fitting procedure as employed for wildtype or control cells, resulting in different parameter sets that can be plugged into the network to assess their implications for network behavior. Likewise, this could be done for the synaptic parameters using paired recordings and recent methods to fit the parameters of short-term synaptic plasticity models to these data [79].
In summary, we have provided a prefrontal cortex network model here with single cells and synapses strictly parametrized through in vitro electrophysiological findings (no specific tuning or adjustment of synaptic currents to compensate for simulated network size), with realistic cellular and synaptic heterogeneity, and with a structural layout derived from anatomical data. We have then systematically compared the full network activity to a number of spiking and correlation statistics from in vivo multiple single cell recordings in awake rodents, as well as LFP data from these animals, and estimates of membrane potential fluctuations from in vivo patch-clamping. Our model is therefore highly validated at the in vivo physiological level, yet it is computationally efficient by virtue of its computationally comparatively simple single unit design. We therefore hope that this network model can serve as a valuable tool in the further study of how physiological and anatomical properties relate to cortical network dynamics, and ultimately cognition, and how alterations of these properties may give rise to symptoms observed in various psychiatric conditions.

Model specification
Neuron model. Single neurons were modeled by the simplified adaptive exponential integrateand-fire neuron (simpAdEx) introduced in [15]: where C is the membrane capacitance, g L a leak conductance (with reversal potential E L ), τ m and τ w are the membrane and adaptation time constants, respectively, Θ denotes the heavyside function, and w V is the V-nullcline of the system as defined in Eq 1. Like the full AdEx [80], this model consists of one differential equation for the membrane potential V (including an exponential term with slope parameter Δ T , which causes a strong upswing of the membrane potential once it exceeds V T ), and one for an adaptation variable w, and can reproduce a whole variety of different spiking patterns [15]. A spike is recorded whenever V crosses V up , at which point the voltage is reset to V r and spike-triggered adaptation is simulated by increasing w by a fixed amount b. The simpAdEx was derived from the full AdEx based on phase-plane considerations, effectively dissecting the dynamics into three different regimes (defined through their distance from the V-nullcline, see Eq 2), each of them approximated in a way that allows for closed-form expressions for the instantaneous and steady-state firing rates. This enables fast and efficient fitting of the model to f-I and I-V curves as commonly used to characterize the electrophysiological behavior of cells in vitro [15] (Fig 1A).
We had shown previously that this model, although estimated from f-I and I-V curves only, can predict spike times under in vivo-like conditions with high accuracy from physiological recordings not used for model fitting [15]. Different from [15], the upper voltage limit V up was initially estimated from the inflection point of the voltage traces. This makes V up an absolute firing threshold (as in the leaky integrate-and-fire neuron) and leaves V th as a free parameter for the subthreshold dynamics, resulting in a shallower exponential rise to the spike, more akin to what would be expected from the action of persistent sodium channels [81] or L-type calcium channels [82].
We estimated neuron models for a large number of in vitro recordings from different cell types from the prefrontal cortex of rats and mice, namely layer 3 (n = 34) and layer 5 pyramidal cells (n = 108), fast-spiking (n = 32), and bitufted (n = 22) interneurons. Additionally, we extracted statistics (means and variances) about f-I curves and subthreshold dynamics of Martinotti cells from the literature [26,[83][84][85][86][87][88], and used these to construct 100 sets of f-I and I-V curves drawn from Gaussian distributions instantiated by the empirically estimated parameters. For each data set drawn from these distributions, Martinotti cell models were estimated. The pool of estimated models for each cell type defines a multivariate parameter distribution for each type of neuron, from which the final parameter sets for the 1000 neurons used in the network simulations were drawn. This joint parameter distribution for each cell type was initially modeled as a multivariate Gaussian, where marginal distributions not of Gaussian shape (as estimated from the empirical data) were first Box-Cox-transformed to adhere with the Gaussian assumptions. In a second step, the Box-Cox transform was inverted to regain the non-Gaussian shape of the marginal distributions (red curves in Fig 1C). The mean values and standard deviations of all model parameters for the different cell types are given in Table 1.
Network anatomy and connectivity. The network is divided into two laminar components, representing the superficial layers L2/3 and the deep layer L5 (Fig 2A). The network also includes a horizontal organization into distinct columns which are typically about 300μm wide [89], so the model is in principle suited to study information transfer between columns. For most part, however, the present analyses focuses on a single column which was found to be sufficient to reproduce in vivo-like resting-state activity, provided a source of constant external input (see below). The relative numbers of pyramidal cells and interneurons in each layer were taken from [58] who studied the rat motor cortex, as such data are not available for the PFC. Following [58], 47% of all cells were modeled as L2/3 pyramidal cells (L2/3-E), 10.4% as L2/3 interneurons (L2/3-I), 38% as L5 pyramidal cells, and 4.6% as L5 interneurons. With regards to the specific types of interneurons and their distribution across layers, we followed [90] and [91] and defined local interneurons (IN-L) with projections within the same layer and column as fast-spiking cells, cross-layer interneurons (IN-CL) as bitufted cells, and far-reaching interneurons (IN-F) with projections both outside of their column and layer of origin as Martinotti cells [90] ( Table 1). The cross-column cells (IN-CC) have been classified as large basket cells [90,92,93], with electrophysiological properties resembling those of pyramidal cells [94]. Therefore, we used the same parameter distributions as for the pyramidal cells in the respective layer for this cell class. Markram et al. [90] also estimated the relative numbers of different types of interneurons within each cortical layer. Together with the classification above, these data result in the full distribution of cell types summarized in Table 2.
Neurons were randomly connected with distinct connection probabilities p con for each pair of cell types as derived from a survey of about 40 studies, e.g. [95][96][97][98][99], most of which are reviewed in [22] and [27], except [18][19][20][21][23][24][25] and [26]. Most of them performed whole cell or dual sharp electrode recordings in vitro in various neocortical regions of rats and mice. We also included a few studies using monkeys, ferrets and cats, as there are more studies from PFC in these species and some parameters were not available in rodents. Connection probabilities were further adjusted jointly with the connection weights to match data from photostimulation experiments as explained in detail below [57,70,71,100,101]. Pyramidal cells within the same layer form clusters of increased connection probability as defined by the "common neighbor rule" [28,29] which states that the connection probability of two neurons increases linearly with the number of neurons they are both connected to. Furthermore, a fraction of 47% of the connections was specified as reciprocal [32], since the proportion of reciprocal connections was experimentally observed to be significantly higher than chance [32,37]. For cross-column projections, connection probabilities exponentially decay with the distance from the column of origin. Data from rodent studies [29,60,77,102,103] suggest a wide range of spatial decay constants. We use the median from these studies, which is 114μm for pyramidal cells and 95μm for interneurons.
Apart from the recurrent synaptic connections within the network, we also introduced constant background currents that are fed into all neurons and which differ in strength for pyramidal cells and interneurons in layer 2/3 and 5. Appropriate values for these four streams of background inputs were determined using reduced equivalent-population input models (see below). It is emphasized that there was no source of external noise fed into the network, i.e. the external inputs consisted of just constant (DC) currents. Thus, all variability observed in the network arises from its internal dynamics.
Synaptic properties. Neurons were connected through conductance-based AMPA-, GABA A -, and NMDA-type synapses, with kinetics modeled by double exponential functions [104] where X 2 {AMPA, GABA A , NMDA}. The reversal potential E rev is set to zero for AMPA and NMDA, and to −70 mV for GABA A [19,24,55,61]. The onset and offset time constants τ on and τ off are set to 1.4 ms and 10 ms, respectively, for AMPA [39,94], 3 ms and 40 ms for GABA A and 4.3 ms [105] and 75 ms [38,39] for NMDA. NMDA conductances exhibit a nonlinear voltage-dependency s(V) due to their magnesium block at lower voltages [106]. Synaptic transmission delays τ D were drawn from Gaussian distributions with means and standard deviations depending on the pair of connected cell types, with parameters derived from the same electrophysiological literature as the connection probabilities (see below; Table 3). Synaptic delays were chosen to increase linearly with distance from the target column [107], τ D (d) = τ D (1 + d), where d is the number of columns separating the connected neurons. Synapses were also equipped with short-term plasticity dynamics implemented by the corrected version [108] of the Tsodyks and Markram model [30] These recursive equations describe the dynamics of the relative efficiency a(t sp k ) across series of spikes, with initial conditions u 1 = U and R 1 = 1, where t sp k is the interval between the (k − 1)th and the kth spike. Model parameters U, τ rec and τ fac were specified according to [31] and [32] who differentiated between facilitating (E1/I1), depressing (E2/I2) or combined (E3/I3) shortterm dynamics, for both excitatory (E) and inhibitory (I) connections (Fig 2B, right panel; Table 4). The cell types of the pre-and postsynaptic neurons determine which of these classes is used for each individual combination (Fig 2B, left panel). Synaptic inputs were further subject to release failures with a probability of 30% [33][34][35][36]. Distributions of peak conductances ("synaptic weights") g max for each cell population were derived in two steps. As the first step, initial estimates were obtained from the anatomical and electrophysiological literature (see above). Generally, peak conductances were adjusted such that they reproduced log-normal distributions of postsynaptic potential (PSP) amplitudes as reported in [37] (means and standard deviations given in Table 3). For excitatory synapses, only the AMPA conductances are specified this way, while NMDA conductances are given by 1.09 times the respective AMPA peak conductance [38,39], with both AMPA and NMDA synapses activating after the same delay. For synaptic connections where peak conductances were not directly available from the surveyed literature, estimates were obtained in one of the following ways: 1) Missing estimates for specific interneuron types were replaced by estimates from other interneuron types where possible. 2) Missing estimates for inhibitory connections within one layer were replaced by those from another layer, rescaled such that they followed the same between-layer ratio as the excitatory inputs. 3) If only means but no standard deviations for the distribution of synaptic parameters were available, we used standard deviations from another layer scaled according to the ratio of the means between layers (missing values of connection probabilities or synaptic delays were estimated in the same way). Finally, for cross-column projections, synaptic weights were assumed to decay with the same exponential course (space constant of 114μm for pyramidal cells and 95μm for interneurons) as taken for the connection probabilities themselves (see above).
In a second step, since by far most of the studies cited above have been performed in sensory areas, data from laser scanning photostimulation (LSPS) [57,70,100,101] and genetically targeted photostimulation [71] studies from motor cortex were used to obtain values closer to PFC. Specifically, all connection probabilities and synaptic weights were scaled such that the total input to each cell type would match the one observed experimentally in these studies. To compute this scaling factor s I i;j , the product p con Á g max of connection probability and synaptic weight was assumed to be proportional to the quantity I LSPS obtained in the experimental studies [70], yielding where |Á| denotes the sum over all matrix elements. p con and g max are then multiplied elementwise by ffiffiffi s I p , such that p con Á g max agrees with I LSPS . The scaled values of all parameters are given in Table 3. The average connection strength (as defined by p con Á g max ) between pyramidal cells and interneurons in the different stripes and layers is coded by the arrow width in Fig 2A. Simulation details. All simulations were performed in customized C code written by the authors. Differential equations were numerically integrated using a 2 nd -order Runge-Kutta method with a maximum time step of 0.05 ms, and all spikes, synaptic, and external events were exactly timed by adjusting the time steps accordingly. More specifically, whenever an incoming spike or a change in external currents occurs within the default time step, the time step is reduced accordingly and all equations are updated at the precise time of that event. Neurons were initialized with V i ð0Þ ¼ E i L and w i (0) = 0 for all i. MATLAB-based routines were used for parameter estimation and network analysis. All software is publicly available at https:// www.zi-mannheim.de/index.php?id=626 and on the freely available repository ModelDB (http://senselab.med.yale.edu/ModelDB/).

Estimation of background currents
The constant background currents I used in the simulations are assumed to replace missing synaptic input from the surrounding network not explicitly simulated. A network of sufficient size should be able to produce these amounts of current inherently (with physiologically realistic synaptic efficacies as used here) and thus self-sustain its in vivo-like activity. To test this idea, we computed the amount of current that is produced by a larger network N that is set up and connected exactly the same way as the actually simulated network n and compared it to the range of background currents that are required for in vivo-like activity. The currents I N were modeled as the time-averaged synaptic currents that are elicited in a single neuron for each cell type (with the averaged cellular and synaptic parameters of all neurons of that type) in response to a bombardment of spikes drawn from a Poisson distribution that mimics a large number of input neurons, reflecting its input connectivity. For a Poisson input spike train with the same firing rate as in the original network, this yields the same synaptic current I n as in the full simulation. Larger networks N are simulated by using a higher number of inputs, resulting in higher overall input spike rates. Because connection probabilities decay with distance between neurons [89], we independently tested two ways to increase the number of input neurons in the network: By increasing its spatial size L (measured in number of columns N C , each L C = 300μm in diameter) or its within-column neuron density D. For cross-column input, we assumed a radial distribution of inputs [60,89] and an exponential decay of connection probabilities with distance (l E = 114μm, l I = 95μm, see above). Only pyramidal cells in L2/3 as well as cross-column (IN-CC) and far-reaching (IN-F) interneurons project across columns. Specifically, the number of input neurons N i syn projecting onto a neuron of cell type i is given by ð8Þ p i con is the fraction of neurons that connects to cell type i. N i D ðL 1 ; L 2 Þ denotes the number of neurons within a hollow cylinder defined by the inner radius L 1 and the outer radius L 2 and the height of 1 mm [60] that are connected to a neuron of cell type i at the center of this cylinder, given a neuron density D. Thus, the two terms represent input from within the same column (up to L C ) and from outside that column (up to the full radius L). The ratio of pyramidal cells and interneurons that project beyond a single column (Table 2) is reflected in different scaling factors for excitatory and inhibitory connectivity for input from within (f C ) and from outside the same column (f N ). The resulting background currents I N are then computed independently for the four main cell types-pyramidal cells and interneurons in layer 2/3 and 5.

Analysis of simulated and in vivo data
Two in vivo data sets were used for comparison with simulation results (kindly provided by Dr. Christopher Lapish, Indiana University Purdue University, Indianapolis and Dr. Thomas Hahn, Central Institute of Mental Health and BCCN Heidelberg-Mannheim). For spike trains and local field potentials, extracellular multiple single-unit recordings were obtained from the rat's anterior cingulate cortex (ACC) while they were performing an eight-arm radial maze task [40]. Stationary periods (largely free from motor or sensory responses) were obtained from 381 units using a previously described stationarity-segmentation method [41]. For the voltage traces, we used patch-clamp recordings from anaesthetized rodents (see [109,110] for details).
Spike trains, voltage traces and local field potentials from the network simulation and the in vivo data were analyzed the same way. For each model cell or recorded unit, mean and coefficient of variation (C V ) of the interspike interval (ISI) distributions were computed. Autocorrelations of ISI series and the zero-lag cross-correlation CC(0) between ISIs from pairs of spike trains were computed as well according to the procedures described in [41] to correct for nonstationarities. All analyses were restricted to spike trains of at least 10 spikes to yield sensible estimates of single-cell statistics without cutting off too much of the low-rate tail from the distributions (see Results section for a more empirical justification).
Similarity among simulated and experimentally obtained distributions was assessed by twosample Kolmogorov-Smirnov (KS) tests, where test statistic D KS is bounded between zero (complete overlap) and one (maximally dissimilar distributions). Underlying distributions were inferred through kernel-density estimation [111] (implemented by the function "ksdensity" in MATLAB's statistics toolbox). As KS test statistics may depend on sample size but different simulations vary in number of spikes, we limited the number of data points to a sufficiently low, common value (30), repeated KS tests with 100 random drawings, and report averages across obtained p values and KS statistics. The overall similarity of a simulation data set with the experimental spike data is quantified by conducting the test for the mean ISI, the C V and the CC(0) distributions, and reporting the minimal p or the maximal D KS value of those three (i.e., the value associated with the largest difference between the compared distributions).
Finally, we visualize the statistical overlap of distributions by plotting shaded areas representing the SEM around the mean at each value of the experimental distributions, which are computed from the 100 bootstrap samples as indicated above.