Network Models Predict That Pyramidal Neuron Hyperexcitability and Synapse Loss in the dlPFC Lead to Age-Related Spatial Working Memory Impairment in Rhesus Monkeys

Behavioral studies have shown spatial working memory impairment with aging in several animal species, including humans. Persistent activity of layer 3 pyramidal dorsolateral prefrontal cortex (dlPFC) neurons during delay periods of working memory tasks is important for encoding memory of the stimulus. In vitro studies have shown that these neurons undergo significant age-related structural and functional changes, but the extent to which these changes affect neural mechanisms underlying spatial working memory is not understood fully. Here, we confirm previous studies showing impairment on the Delayed Recognition Span Task in the spatial condition (DRSTsp), and increased in vitro action potential firing rates (hyperexcitability), across the adult life span of the rhesus monkey. We use a bump attractor model to predict how empirically observed changes in the aging dlPFC affect performance on the Delayed Response Task (DRT), and introduce a model of memory retention in the DRSTsp. Persistent activity—and, in turn, cognitive performance—in both models was affected much more by hyperexcitability of pyramidal neurons than by a loss of synapses. Our DRT simulations predict that additional changes to the network, such as increased firing of inhibitory interneurons, are needed to account for lower firing rates during the DRT with aging reported in vivo. Synaptic facilitation was an essential feature of the DRSTsp model, but it did not compensate fully for the effects of the other age-related changes on DRT performance. Modeling pyramidal neuron hyperexcitability and synapse loss simultaneously led to a partial recovery of function in both tasks, with the simulated level of DRSTsp impairment similar to that observed in aging monkeys. This modeling work integrates empirical data across multiple scales, from synapse counts to cognitive testing, to further our understanding of aging in non-human primates.


INTRODUCTION
Spatial working memory in rhesus monkeys, as in humans, is mediated by the action potential firing activity of neurons in the dorsolateral prefrontal cortex (dlPFC) (reviews: Funahashi, 2017;Miller et al., 2018). And, as in humans (Albert, 1993;Salthouse et al., 2003;Fisk and Sharp, 2004;Rhodes, 2004;Sorel and Pennequin, 2008), the effects of aging on working memory is heterogeneous-while a significant proportion of rhesus monkeys become increasingly impaired on tasks of spatial working memory during normal aging ("unsuccessful agers"), a significant number of "successful agers" show no signs of impairment (Lacreuse et al., 2005;Moore et al., 2006Moore et al., , 2017. Furthermore, for impaired subjects, impairment typically begins fairly early in the aging process, during early middle age (Moore et al., , 2006(Moore et al., , 2017. This distribution of impaired and spared monkeys across the lifespan enables assessment of brain changes associated not simply with aging, but with cognitive performance per se in this model of normal human aging. This approach has been taken in studies of normal aging by our group and others (reviews: Peters, 2007;Luebke et al., 2010Luebke et al., , 2015Wang et al., 2011;Peters and Kemper, 2012), in which monkeys were assessed on working memory tasks such as the Delayed Response Task (DRT) and the Delayed Recognition Span Task in the spatial condition (DRSTsp), and their brains subsequently examined for a wide variety of parameters. Thus, declines in spatial working memory have been associated with many sub-lethal changes to the structure and function of neurons, glial cells, and white matter pathways in the dlPFC of the rhesus monkey (Peters, 2007(Peters, , 2009Peters et al., 2008;Bowley et al., 2010;Luebke et al., 2010;Shobin et al., 2017). Reductions in synapses and increased dystrophy of white matter pathways begin in early middle age; for example, Peters et al. (2008) showed a continuous decrease in the number of excitatory and inhibitory synapses, detectable even in middle age in the monkey dlPFC.
In addition to well-documented structural changes, functional alterations to the electrophysiological properties of supragranular neurons in the aging monkey dlPFC have been reported in several studies (Chang et al., 2005;Wang et al., 2011;Coskren et al., 2015). Layer 3 dlPFC pyramidal neurons exhibit significantly increased evoked action potential (AP) firing rates in vitro for aged compared to young monkeys (Chang et al., 2005;Coskren et al., 2015). In contrast, another study (Wang et al., 2011) reported that in vivo firing rates of dlPFC DELAY neurons decreased with aging during the DRT, for monkeys performing the task well. Since persistent firing patterns of dlPFC pyramidal neurons encode precisely tuned spatial and temporal information during a working memory task, exploring how these findings might complement each other may reveal new predictions about spatial working memory decline.
While the growing body of literature over the past 20 years has documented numerous changes to neurons, glial cells and white matter pathways in the aging brain that are associated with cognitive changes, which changes are the key determinants of age-related cognitive declines has not yet been firmly established (Hof and Morrison, 2004;Luebke et al., 2010;Morrison and Baxter, 2012;Peters and Kemper, 2012;Konar et al., 2016;Motley et al., 2018;Cleeland et al., 2019). This lack of insight is partially due to the difficulty in first, selectively targeting individual variables (such as firing rate or synapse number) in humans or experimental subjects without, second, also altering upstream or downstream effectors. A powerful way to gain insight into which age-related changes are most consequential for dlPFC network behavior during normal aging is through computational models that are constrained by empirical data. This is the approach we use here through a systematic study of the "bump attractor" neural network model, and its extension to memory retention in a new task similar to the DRSTsp.
The bump attractor network has been used to model spatial working memory in the DRT (Compte et al., 2000;Wang et al., 2011Wang et al., , 2013Wimmer et al., 2014;Wu et al., 2016 for review). The bump attractor exhibits persistent firing activity after an initial stimulus disappears, encoding a memory of the stimulus location. The bump attractor was used to show that an age-related loss of synaptic strength could account for reduced firing rates of dlPFC DELAY neurons (Wang et al., 2011), but the effects of increased excitability of pyramidal neurons seen in vitro on the function of bump attractor models has not been examined. We have previously used computational modeling to predict which intrinsic electrophysiological and morphological properties of individual pyramidal neurons contribute to the action potential firing rate increases seen in aging (Coskren et al., 2015;Rumbell et al., 2016). Here, we extend these previous single neuron studies using network-level modeling to predict the effect of empirically observed physiological and structural changes in the aging rhesus dlPFC on cognitive behavior. We obtained whole-cell patch clamp recordings of dlPFC pyramidal neurons from behaviorally characterized monkeys across the adult lifespan. Empirical results indicated that previously reported physiological changes seen in aging are already present in middle age, and are correlated with cognitive impairment. This modeling work aims to connect age-related changes in non-human primates across multiple empirical scales, from synapse counts and physiology of single neurons up to network output and cognitive performance. We predict that the increased firing rates and reduced synapse density observed with aging may partially compensate one another, but still are sufficient to induce substantial cognitive deficits in the DRT and DRSTsp.

Experimental Subjects
The rhesus monkeys (Macaca mulatta) studied here were a part of a larger study of normal aging of the brain. Data were obtained from a total of 9 young (8.9 ± 0.5 years old), 18 middle-aged (18.2 ± 0.5 years old), and 12 aged (24.7 ± 0.7 years old) monkeys ( Table 1). Prior to the electrophysiological experiments, all but 3 young monkeys completed cognitive testing of the DRSTsp (Moore et al., 2017). Monkeys were obtained from the Yerkes National Primate Research Center at Emory University (Atlanta, GA, USA), and housed in the Laboratory Animal Science Center (LASC) at the Boston University School of Medicine (BUSM) under strict accordance with the guidelines established by the

Behavioral Assessment
A total of 36 subjects were tested on the DRSTsp, 6 young, 18 middle-aged, and 12 aged monkeys. The DRSTsp is a short-term memory task in which the subject distinguishes a novel cue (spatial location) from an increasing set of recently presented, familiar cues (spatial locations) (Moore et al., 2017 and Figure 1A). All stimuli were identical brown disks, each covering one of 18 small wells arranged on a testing board in a 3 × 6 matrix. On the first trial one of the wells was baited with a reward then covered with a disk. Once a screen was raised, the monkey moved the disk to receive the reward. The screen was lowered before the second trial, when that same well (with no reward) was covered with a disk and a second well was baited with a reward and covered with an identical disk (10 s delay). Subsequent trials include the previous cue plus one new cue, until the monkey is unable to identify novel cues (spatial locations). Monkeys were tested for 10 consecutive days. The cognitive score recorded is the recognition span: the mean number of unique cues the monkey identified correctly before making a mistake. Some of the subjects were assessed on the DRSTsp more than once over a period of several years as part of a longitudinal study of aging. The DRSTsp score used in the present study for longitudinally assessed subjects was the one most proximal to the date the subject was sacrificed.
Whole-cell patch clamp recordings from 324 layer 3 dlPFC pyramidal neurons from 38 subjects were used in this study. In total, 40 neurons for the young, 188 for the middle-aged and 96 for the aged subjects.

The Bump Attractor Network Model
The "bump attractor" network model was originally designed to simulate a spatial working memory task for awake behaving monkeys known as the DRT or the Oculomotor Delay Task (Figure 2A, Goldman-Rakic, 1995;Compte et al., 2000;Wang et al., 2011;Wimmer et al., 2014). After the monkey fixates on the center of a computer screen (for 1 s in our simulations), a light appears during the cue period in one of eight directions around a circle (occurring at 0 • and lasting 0.5 s in our simulations). After the cue period the light turns off, and the monkey fixates on the center of the screen during a delay period (2 s in our simulations). During the response period (0.7 s in our simulations), the monkey shifts its gaze (makes a saccade) to indicate its memory of the location of the stimulus. The selective firing of pyramidal neurons in the dlPFC to different spatial locations (their preferred directions) is thought to be the mechanism by which stimulus location is encoded. Thus, monkeys performing this task successfully have neurons that exhibit persistent neural activity during the delay period tuned to the stimulus location (Goldman-Rakic, 1995;Compte et al., 2000;Wimmer et al., 2014). This is achieved in the bump attractor model with excitatory neurons tuned to each spatial location (so that we can label each neuron by that location-its "preferred direction"), with the synaptic strength between pairs of excitatory neurons determined by the difference between the neurons' preferred directions.
The bump attractor network used here assumed a firingrate model for individual neuronal dynamics, modified from Wimmer et al. (2014). The code was written in Matlab and is available on ModelDB (McDougal et al., 2017) at Accession #256610. The local cortical network was composed of N E = 640 excitatory (pyramidal) neurons (80%) and N I = 160 inhibitory neurons (interneurons, 20%) (Abeles, 1991;Braitenberg and Schütz, 1991) (Figure 2B). Each excitatory and inhibitory neuron received three general types of synaptic inputs: from excitatory and inhibitory connections with other neurons in the same network, and excitatory connections from other unspecified areas external to the network. The modeled excitatory neurons also received an external input current representative of the "cue." The synaptic excitation of each neuron was modeled with two distinct equations, separately representing synaptic inputs mediated by AMPA and NMDA receptors. The voltage dependence of the NMDA currents was omitted from the present study, since Compte et al. (2000) showed that the slow kinetics of NMDA-mediated synaptic transmission was the most important feature for generating stable persistent activity.
The equations for FR of each excitatory and inhibitory neuron (r E and r I ) were given by for j = 1, . . . , N E and k = 1, . . . , N I . The parameters τ E and τ I represented the membrane time constants for excitatory and inhibitory neurons, respectively, (in ms), with f (I) the FR activation function (or f-I curve, in Hz), total synaptic input currents I E and I I (in pA), and Gaussian white noise inputs σ E ξ j E (t) and σ I ξ k I (t) (in Hz), where σ E and σ I were the standard deviations. In this study, τ E = 20 ms, τ I = 10 ms, σ E = 1, and σ I = 3. The activation function for the firing rate of each neuron took the piecewise form (Brunel, 2003;Wimmer et al., 2014;Roxin and Compte, 2016): The parameter I c represented a current near rheobase (measured in pA), with corresponding firing rate ν c (in Hz) (see Figure 2C). Thus, the parameter pairs {ν ce , I ce } and {ν ci , I ci } uniquely defined the f-I curves for excitatory and inhibitory neurons.
Here, I ci = 20 pA and ν ci = 50 Hz, values approximately in the regime for Chandelier cells in the monkey PFC (Zaitsev et al., 2009;Povysheva et al., 2013). The parameters ν ce and I ce were chosen to fit our empirical data as described below. See also the Supplementary Information for a general, less excitable form of both f -I curves.
The total synaptic input currents (in pA) for each neuron were with time-dependent synaptic currents given by The subindices a and n represented the AMPA and NMDA synaptic receptors for the glutamatergic synapses, and g the Frontiers in Computational Neuroscience | www.frontiersin.org GABA A receptor. Parameters τ a , τ n , and τ g were the synaptic decay time constants for AMPA, NMDA, and GABA A (Hestrin et al., 1990;Spruston et al., 1995;Salin and Prince, 1996;Xiang et al., 1998). We set τ a = 2 ms, τ n = 100 ms, τ g = 10 ms. The quantity r represented the mean FR value of each neuron type, and the G pq parameters with {p, q} = {E, I} were the constant synaptic weights. The connectivity matrix M was a circular Gaussian function describing the translation-invariant connections among the excitatory neurons. Labeling each neuron by its preferred direction, the matrix element M ji for two neurons j and i depended on the difference between their preferred directions. That is, being the preferred stimulus direction of the excitatory neuron j, where a complete circumference was divided into N E angles (see Figure 2D). The quantities I 0E and I 0I represented constant excitatory input currents from other brain areas to excitatory and inhibitory neurons, respectively, and I j stim the input current to each excitatory neuron due to stimuli during cue periods of the tasks described below. Here, I 0E = 80 pA and I 0I = 15 pA. The parameters G pq were varied using the Latin Hypercube Sampling design, as described below.
The input stimulus currents for each excitatory neuron representing a stimulus at 0 • was modeled as in Wimmer et al. (2014) as I j stim = I st e c · cos θ j N E j=1 e c · cos θ j , which provided the bump-shaped input current to all excitatory neurons ( Figure 2D). Parameter I st defined the strength of the stimulus and c the concentration of excitatory neuron to excitatory neuron connectivity. In our simulations, I st = 40,000 pA, and c = 1.5 for the DRT model and c = 20 for the DRSTsp model. A typical simulation of the DRT is shown in Figure 2E.

Increased AP Firing Rate (Hyperexcitability) of Pyramidal Neurons With Aging
We fit parameters ν ce and I ce of the excitatory neurons f-I curve to in vitro FR data from the subjects described above. We computed the mean FR vs. different input currents of all neurons from each monkey, then computed a grand mean of FR vs. input current for monkeys in each age group ( Figure 3B). Fits of the parameters ν ce and I ce to these data were estimated manually using the FindFit function from Mathematica (Wolfram Research, Champaign, IL). The value I ce = 98 pA fit all three age groups well, and also agreed with rheobase values for regular spiking pyramidal neurons with low input resistance (Zaitsev et al., 2012). We used ν ce = 5, 7, and 9 Hz to fit the f-I curve to data from the young, middle-aged, and aged groups, respectively ( Figure 3B).
Thus, to fit in vitro FR data it was sufficient to increase the single parameter ν ce . Assuming the magnitude of changes to f-I curves with aging would be similar at physiological vs. room temperature, we used these parameter values directly in our "young, " "middle-aged, " and "aged" network models. We also repeated some DRT simulations after fitting a generalized, less excitable firing rate activation function to the young, middleaged, and aged monkey data (Supplementary Information and Supplementary Figure 1a). Peters et al. (2008) reported a loss of the number of asymmetric and symmetric synapses (excitatory and inhibitory, respectively) in Layers 2/3 of area 46 of rhesus monkey as adult monkeys age: 10% in middle-aged monkeys, and 30% in aged. For the "young monkey model cohort" DRT studies below we represented these respective losses as a 10% decrease in the synaptic weights G EEa , G EEn , and G IE in middle-aged models, and 30% decrease in these parameters in aged models. For the "young monkey model cohort" DRSTsp studies, we decreased these parameters in a semi-continuous way for the simulated "aged" networks.

Exploring the DRT Model Parameter Space
We conducted two general types of parameter exploration studies for the simulated DRT. The first were sweeps across parameter space using a space-filling Latin Hypercube Sampling (LHS) design, as in (Johnson et al., 1990;Rumbell et al., 2016). The LHS design identified 4,200 points (networks) across the parameter space of synaptic weights (G EEa , G EEn , G IE , G EIa , G EIn , and G II ), representing semi-random combinations of the synaptic weights homogeneously distributed across the 6-dimensional space (Rumbell et al., 2016) (see Table 2). Bounds of the parameter space were chosen so that the mean firing rate of excitatory neurons in all "young" model networks performing the DRT successfully had mean firing rates of excitatory neurons near 35 Hz, a realistic value for the PFC (see Figure 5, black curves). We examined the results of network simulations throughout the LHS as the excitatory firing rate parameter ν ce varied from 2 to 10 Hz (Figure 3). All other parameters remained fixed for all simulations. For each of the 4,200 parameter combinations in the LHS, the DRT was simulated 7 times to account for variations due to noise in the firing rate equations (Equation 1), and results were averaged across all repetitions. We identified five main categories of network behavior: maintaining tuned persistent activity (TPA) to any angle orientation until the end of the delay period; maintaining TPA tuned to the original stimulus location at 0 • (TPA-S) until the end of the delay period; underexcited networks (either unable to maintain TPA until the end of the delay, or unable to generate a bump of activity at all); over-excited networks (all neurons firing at a similar rate by the beginning of the cue period and for the remainder of the task); and partially over-excited networks (initially maintaining TPA, but becoming over-excited sometime during or after the cue period).
For the second type of parameter exploration, we first defined a young model cohort as points in the LHS associated with ν ce = 5 Hz which maintained TPA-S until the end of the delay period in at least one of the 7 simulation repetitions (shown as the ν ce = 5 Hz bar in Figure 3D). We then perturbed parameters of the young model cohort to simulate the effects of morphological and physiological changes that have been observed with aging, or that we propose as novel possibilities (Figure 4). Eight kinds of perturbations were made to the young model cohort: (1) The increased excitability condition: increasing ν ce from 5 Hz to 7 and 9 Hz, keeping all else constant, to simulate the effect of increasing excitability of individual neurons as monkeys age beyond young adult.
(2) The synapse loss condition: reducing the values of G EEa , G EEn , and G IE by either 10 or 30%, simulating the effect of fewer excitatory and inhibitory synapses onto pyramidal neurons for middle-aged and aged monkeys, respectively. (3) Combined increased excitability and synapse loss condition: simultaneously perturbing the excitability and synapse parameters from (1) and (2). Below this perturbation is also called the "observed data conditions." (4) The observed data conditions (combined increased excitability and synapse loss conditions) as in (3), plus increasing the excitability of the inhibitory neurons in the same proportion as in the excitatory neurons: increasing ν ci from 50 Hz to 70 and 90 Hz for middle-aged and aged simulations, respectively. Each of these four perturbations was then repeated after adding short-term synaptic facilitation (as described below) to the networks of the young model cohort. As with the LHS, each young cohort network simulation was repeated 7 times with different randomized seeds and results were averaged.
During in vivo recordings, Wang et al. (2011) found a similar reduction of the firing activity of the DELAY neurons during the cue and delay periods of a DRT with aging, in monkeys performing the task well, and they used computational modeling of the DRT to show that an age-related loss of synaptic strength could account for that reduction. Here we studied if this was still true when the excitability of the excitatory neurons increased at the same time ( Figure 5A). To do this, we examined the mean firing rate of excitatory neurons for the networks in the LHS maintaining TPA-S as ν ce increased from young to aged models (for ν ce = 5, 7, and 9 Hz). In these simulations we reduced the stimulus scaling factor, I st , by 10 and 30% for middle aged and aged networks, respectively. This made firing rates during the cue period more consistent with data from Wang et al. (2011), without significantly affecting firing rates during the delay period (driven primarily by G EEa and G EEn ).

Simulating Memory Retention in a DRSTsp-Like Task
Most of the computational models developed to understand the neural mechanisms underlying working memory are based on delayed memory tasks with just one stimulus to remember during the delay period. A small but growing number of computational studies simulate behavioral tasks that involve the maintenance of several stimuli in memory. Developing models able to maintain the memory of multiple stimuli is challenging, due to interference produced by the simultaneous activation of the many neurons encoding the different stimuli. One way to avoid this problem is to have each stimulus encoded by sparse bumps of activity that does not overlap with firing activity encoding other stimuli (Amit et al., 2003). The levels of inhibition (Edin et al., 2009) and excitation (Wei et al., 2012) in the network are two factors that control the working memory capacity for maintaining multiple stimuli. The balance between them influences the number of stimuli that can be remembered and how bumps fail (fading out vs. merging). Continuous attractor models require fine-tuning of the parameters to control the balance of excitation and inhibition. Such fine tuning can be relaxed by introducing short-term synaptic plasticity. Mongillo et al. (2008) proposed that calciummediated synaptic facilitation could be the neural mechanism underlying working memory. Rolls et al. (2013) showed that synaptic facilitation can make the system more robust to model parameters. Both found that synaptic facilitation increases the multi-item working memory capacity in a discrete attractor network. Mi et al. (2017) also found that working memory capacity increased with the time constant of synaptic depression. Some bump attractor models of the DRT have included synaptic facilitation (Itskov et al., 2011;Hansel and Mato, 2013), but in the past it has often been omitted.
We extended the simulation setup of the DRT to create a model of memory retention in a simplified task that is similar to the DRSTsp. Our original bump attractor model had great difficulty maintaining the memory of multiple stimuli simultaneously, due to interference between the stimulus responses (Amit et al., 2003;Barak and Tsodyks, 2014) and because during the first cue and delay periods, all neurons not involved in the bump activity were highly inhibited. Thus, when introducing subsequent stimuli, the stimulus input current (equal for all stimuli) was insufficient to induce subsequent bumps of firing activity. This difficulty was addressed by adding shortterm synaptic facilitation in the excitatory-to-excitatory synaptic connections. As in Itskov et al. (2011), we assumed that the overall behavior of the synapses was facilitating, and modified the first two equations in Equation (3) as with the facilitation variable u given by Here u is a utilization variable defining the fraction of total neuronal resources that each spike is able to access, representing increased calcium levels at presynaptic terminals (Mongillo et al., 2008). Without firing, u decays to its baseline value U with time constant τ f . In all our simulations U = 0.001 and τ f = 1, 500 ms. Behavioral data showed a maximum DRSTsp score of 3.61 for one of the middle-aged subjects ( Figure 1B). The cognitive score recorded is the recognition span: the mean number of stimuli the monkey identified correctly before making a mistake. Thus, we designed our computational task so that perfect performance gave a simulated DRSTsp score (hereafter called DRST sim ) of four, corresponding to successful encoding of four stimuli equally spaced around the ring from the DRT (Figure 6A). Our DRSTsp model simulated a "pre-cue" period of 1 s, followed by alternating cue and delay periods (lasting 0.5 and 2 s, respectively). The precue and delay periods represented times in the real DRSTsp when the screen was down while appropriate wells were baited. In each subsequent cue period, a new stimulus was presented in addition to any stimuli that were presented before. In the first cue period a stimulus was presented at 0 • ; in the second, stimuli were presented at 0 • and 90 • . In the third cue period stimuli were presented at 0 • , 90 • , and −90 • , and in the fourth stimuli were presented at 0 • , 90 • , −90 • , and 180 • . Stimuli were identical each time at each presentation, with I st = 40,000 pA and c = 20. After the first cue period, if a network simulation encoded all stimuli during cue period n and maintained all n − 1 stimuli throughout the previous delay period (by maintaining a bump for each stimulus with maximum FR of at least 5 Hz), we assumed that the novel stimulus was identified correctly. In this way, our task models the ability of our network to retain the presented stimuli, but does not model the process of choosing which stimulus among all presented is the novel one. Our task also assumes that each trial period lasts the same amount of time, whereas monkeys performing the DRSTsp terminate each successful trial by moving the disk corresponding to their choice before the maximum trial time expires. Sample networks with DRST sim = 4, 3, 2, and 1 are shown in Figure 6B.
With this new DRST sim established, we again used LHS to explore how task performance varied across the parameter space of the six synaptic weights (G EEa , G EEn , G IE , G EIa , G EIn , and G II ; see Table 3) and for ν ce = 5, 7, and 9 Hz. As for the DRT, each LHS had 4,200 points ( Figure 6C); we chose the bounds of the parameter space so that all 4,200 networks encoded at least the first stimulus for all three ν ce values. Simulations were performed only once, since synaptic facilitation made the network performance much less sensitive to random fluctuations.
We also defined a young monkey model cohort for the DRSTsp model to be a group of 40 randomly chosen points in the LHS with ν ce = 5 Hz ( Figure 6D, young case: black bars), half with DRST sim = 4 and the other half with DRST sim = 3. We then "aged" the networks in this cohort by applying increased excitability and synaptic loss conditions, and increasing the excitability of the inhibitory neurons, as follows. We set bounds for ν ce just beyond our young and aged values (4.5 and 9.5 Hz) and generated 40 randomly distributed values between them. We then defined two piecewise linear functions that interpolated the synapse loss conditions from the DRT (0, 10, and 30% loss when ν ce = 5, 7, 9 Hz, respectively, slope 5% loss/Hz when ν ce < 7 Hz and 10% loss/Hz when ν ce ≥ 7), and computed a corresponding level of synapse loss for each of the 40 ν ce values. We computed corresponding values for ν ci by multiplying each ν ce value by 10, consistent with the relationship between these two parameters elsewhere in our study. We computed an associated age for each ν ce value similarly, interpolating the mean age for the three empirical groups (8.9, 18.2, and 24.7 years for young, middle aged, and aged groups) when ν ce = 5, 7, 9 Hz, respectively. We induced the increased excitability and synapse loss conditions of these 40 values separately and together, plus increased excitability of inhibitory neurons, and then computed the corresponding DRST sim for each transformed member of the young cohort (6,400 simulations in all: 160 for each point in the cohort). We then calculated the mean DRST sim for each simulated age group (young group, 7.5 < years < 13; middle-aged group, 13 ≤ years < 21; and aged group, 21 < years < 26). To summarize the change in performance across the simulated age span (Figure 7), we calculated the relative change in the mean DRST sim from the mean ages of the young simulations to the aged simulations (9.9 and 24.2 years, respectively) under each perturbation condition.

Statistics
All analyses were performed in RStudio 1.1.463 (RStudio Team, 2015;Mangiafico, 2016). We treated age as a continuous variable in most analyses, using standard linear regression as well as generalized linear mixed-effects models (GLMMs) to explore relationships between output variables vs. dependent variables. Because electrophysiological data were collected from several neurons per subject, the subject was treated as a random effect blocking factor in relevant analyses (Grafen and Hails, 2002;Darian-Smith et al., 2013). For the GLMMs, marginal R 2 captured the variance explained by fixed factors; conditional R 2 captured the variance explained by the whole model including the random effects (Nakagawa and Schielzeth, 2013). In analyses with DRSTsp score (one data value per subject) vs. physiological variables, physiological variables were first averaged within each subject. Age was treated as a categorical variable (young, middleaged, and aged) when analyzing relationships firing rate vs. age ( Figure 1E). Specifically, we used two-way repeated measures ANOVA (including subject as a random effect) to determine whether there were differences in FR vs. age group, with injected current as the repeated measure. Before performing the ANOVA, quantile-quantile plots were used to identify extreme outliers from normality for each age group and level of injected current, and then Shapiro-Wilk's test of normality was applied. Post-hoc analyses were conducted using Tukey-adjusted comparison of estimated marginal means. For all analyses, the significance level was α = 0.05.

Rhesus Monkeys Show DRSTsp Impairment and Increased AP Firing Rates of dlPFC Pyramidal Neurons With Aging
We examined relationships between performance on the DRSTsp, age of the monkeys, and physiological variables assessed (Figure 1). Among the 39 subjects analyzed for this study, 36 of them completed behavioral testing on the DRSTsp (6 young, 18 middle-aged, and 12 aged). There was a significant decrease in performance on the DRSTsp with age [linear regression, DRST = −0.032 * Age + 2.97, R 2 = 0.18, F (1, 34) = 7.43, p = 0.0101; Figure 1B] replicating past findings on impaired working memory with age in monkeys (Herndon et al., 1997;Moore et al., 2006Moore et al., , 2017Wang et al., 2011). In vitro slices were prepared from the dlPFC of the same monkeys that were tested behaviorally for whole cell patch clamp recordings. Data from 324 neurons met our criteria for inclusion. Our statistical analysis used random effects modeling to control for the variation of physiological variables within monkeys (Materials and Methods). Input resistance (R n ) increased with aging (R n = 4.02 * Age + 89.1 M , marginal R 2 = 0.08, conditional R 2 = 0.20, p < 0.0001; Figure 1C), while resting membrane potential did not (p = 0.11). There was a significant positive relationship between input resistance and AP firing rate at all levels of current injection (p < 0.00016 and conditional R 2 ≥ 0.21 for each level, Figure 1D).
Next, we compared depolarizing current step-evoked AP firing rate in neurons from young, middle-aged, and aged subjects ( Figure 1E). Before using the two-way repeated measures ANOVA, outliers from normality were removed for AP  In all graphs, young, middle-aged, and aged subjects are shown in black, blue, and red, respectively. Data were obtained from a total of 9 young, 18 middle-aged, and 12 aged monkeys ( Table 1). All but 3 young monkeys completed cognitive testing of the Delayed Recognition Span Task in the spatial condition. Whole-cell patch clamp recordings from 324 pyramidal neurons from 38 subjects were used.
firing rate vs. current injection level for each age group. There was a significant effect of age group on firing rate [F (2, 22) = 8.63, p = 0.0017] but no significant interaction between age group and injected current [F (2, 799) = 2.76, p = 0.064]. Firing rates were significantly higher for the aged vs. young monkeys [Tukeyadjusted comparison of estimated marginal means, t (22) = 4.14, p = 0.0012]; firing rates for middle-aged monkeys did not differ from the other two groups [young vs. middle-aged, t (22) = 2.23, p = 0.072; middle-aged vs. aged, t (22) = 1.99, p = 0.139; Figure 1E]. Hereafter, we refer to this increased firing rate of current stepevoked APs as hyperexcitability of pyramidal neurons.
Given significant relationships between aging and both DRSTsp performance and the physiological variables, we next examined whether both aging and physiological variables might together affect DRSTsp performance. Since R n and firing rate were significantly related at each injection level, we only included a single physiological variable in the analysis: firing rate in response to the 130 pA current injection. There was a significant negative linear relationship between DRSTsp span and firing rate [DRSTsp = −0.055 * FR + 2.95 Hz, R 2 = 0.19, F (1, 33) = 7.71, p = 0.0090, Figure 1F]. To examine the effect of age and firing rate on DRSTsp span, we fitted a regression model that included age, firing rate, and their interaction (R 2 = 0.26, F (3, 31) = 3.68, p = 0.0225). The negative effect of the firing rate on DRSTsp was significant (p = 0.0441) and it did not depend on the age of the animal (age by firing rate interaction effect p = 0.1340). The effect of age alone was not significant either (p = 0.0904). These analyses lend support to the hypothesis that firing rate affects DRSTsp performance more than age does, though this support is somewhat tenuous for the following reasons. First, the delay between DRSTsp assessment and physiological recording of neurons from these valuable subjects can last several months to a year or more, and may mask a relationship between these variables. Second, some of these subjects were employed in longitudinal studies for other projects and tested on the DRSTsp a few times; we cannot rule out a practice effect in these subjects. Thus, we turned next to computational modeling to create a controlled environment for exploring how changes in AP firing rates (independent of a subject's age) might affect performance on spatial working memory tasks.

An Optimal Level of Pyramidal Neuronal Excitability Maximized Successful Performance in the DRT Model
Here, we used computational modeling to explore how these empirical results of physiological changes and a past study on reduced synapse counts with aging might together impact cognitive function in middle-aged and aged monkeys. We chose a continuous bump attractor network model which has been used by several previous research groups (Compte et al., 2000;Wang et al., 2011Wang et al., , 2013Wimmer et al., 2014). The network model simulates an empirically measured behavior (the oculomotor DRT administered to macaque monkeys) through a simple but elegant model network that can be manipulated in specific ways (Figure 2; Materials and Methods). During the cue period of the DRT simulations, all excitatory neurons received an external input current, tuned to 0 • direction, proportional to their orientation preference ( Figure 2D). In the subsequent delay period, the main input to model neurons came through recurrent activity within the network. If persistent neural activity is maintained at the end of the delay period-in the form of a bump-shaped FR profile for a subset of the excitatory neuronsthen the orientation angle corresponding to the center of the bump of activity at the end of the delay is assumed to represent the stimulus location encoded by the network (Figure 2E).
Varying model parameters created a distinct network configuration, and each simulation represented a model DRT trial. Successful DRT trials maintained TPA (and in particular TPA-S) until the end of the delay period (Figures 2E, 3A2). Failed DRT trials can occur in several ways, including under-excitation, a loss of TPA during the cue or delay periods ( Figure 3A1); generating TPA early in the cue or delay periods but ending with full network over-excitation ( Figure 3A3); and an over-excited network throughout the simulation ( Figure 3A4).
An important feature of the bump attractor model we chose to use (Wimmer et al., 2014) is the ability to control the f-I curve of individual neurons by two simple parameters. Adjusting the parameter ν ce was sufficient to fit the f -I curves of excitatory model neurons to the mean in vitro AP firing rates of pyramidal neurons reported above ( Figure 3B and Equation 2). Young, middle aged, and aged f -I curves used ν ce values of 5, 7, and 9 Hz respectively. These fits were accurate for low input current, but became more excitable than the empirical f-I curves for high inputs. However, our results did not depend on the specific shape of the f-I curves: similar results were found when using less excitable analogs of the firing rate function (Supplementary Figure 1). Synaptic weight parameters represented the number of excitatory and inhibitory synapses onto excitatory (G EEa , G EEn , and G IE ) and inhibitory (G EIa , G EIn , and G II ) neurons. Varying these synaptic weights allowed us to model published data showing an age-related reduction in the excitatory and inhibitory synapses with aging in the neuropil of layer 2/3 of rhesus monkey dlPFC (Peters et al., 2008). In a series of simulations, we kept ν ce fixed at distinct values and then used the Latin Hypercube Sampling (Materials and Methods) to choose points throughout the parameter space generated by the six synaptic weights. Each point was simulated seven times, as if administering multiple DRT trials to each virtual monkey. We examined whether points across the parameter space maintained TPA throughout the delay period, or had a different outcome ( Figure 3C).
Increasing the excitability of model excitatory neurons strongly affected the ability of the DRT model networks to maintain TPA. The 3D plots in Figure 3C summarize DRT model results for 4,200 points across the parameter space of synaptic weights as ν ce increased. Blue circles represent "underexcited" networks, where no TPA was detected at the end of the delay period (e.g., Figure 3A1). Green circles represent networks that maintained TPA for the entire delay period (e.g., Figure 3A2)-corresponding to "virtual monkeys" that encoded some stimulus location. Remaining circles represent partially over-excited networks (yellow; e.g., Figure 3A3), which started Ring connectivity among excitatory neurons. M ji is the connectivity matrix element between excitatory neurons j and i, depending on the difference between their preferred directions (θ j and θ i , respectively). Each neuron is strongly connected to its nearest neighbors (red line), with the connection strength decaying with distance (yellow and green lines). The neuron with preferred direction 0 • receives the strongest stimulus current, which follows a Gaussian distribution centered at 0 • , generating a bump of activity during the delay period (persistent activity tuned to the stimulus location) when excitation and inhibition in the network are well-balanced. (E) Sample DRT model output (side and top views) when the 0 • stimulus location is encoded correctly: excitatory neuron FR vs. simulation time with neurons labeled by their preferred direction. out with TPA that led to firing in all neurons, or were already over-excited networks during the fixation period before the cue ever appeared (red; e.g., Figure 3A4). Below we focus on the green points, which maintained TPA throughout the delay period.
The greatest likelihood of maintaining TPA, and in particular TPA-S, across the synaptic weight parameter space occurred for intermediate values of ν ce -neither too low nor too high ( Figure 3D). The optimal value among those we tested was ν ce = 5 Hz. Values of ν ce below this optimum led to many networks across parameter space that were under-excited; above the optimum many networks were over-excited ( Figure 3C). In either case, the networks were much less likely to encode the stimulus successfully until the end of the delay period, and this reduced likelihood was visible even in our "middleaged" networks (ν ce = 7 Hz). Thus, if the excitability of single pyramidal neurons increased above the optimal value without complementary changes in other parameters, the network's ability to maintain TPA, in particular to the stimulus location-and successful working memory performance-was impaired. These findings also held with the less excitable form of the f -I curves fit to the data (Supplementary Figures 1b,d).

Synaptic Changes Partially Compensated for Increased Pyramidal Neuron Excitability in Successful DRT Model Networks
Points across parameter space that did maintain TPA-S as ν ce increased had lower values of the excitatory synapse weights for excitatory neurons (G EEa and G EEn ) and higher values of the inhibitory synapse weights (G IE ) for these neurons (see Figure 3F; standard error bars lie beneath the symbols). Figure 3E shows the mean values of each synaptic weight for all sampled points maintaining TPA-S as ν ce increased. There was no clear relationship between ν ce and the synaptic weight parameters onto interneurons (G EIa , G EIn , and G II ). However, excitatory parameters G EEa and G EEn decreased noticeably as ν ce increased before leveling off after the middle-aged ν ce value, while G IE increased in a nearly linear fashion. To maintain TPA in the (2) TPA-S; (3) partial over-excitation; and (4) over-excitation. (B) The model f-I curve fit (solid lines) to empirical AP firing rates of pyramidal neurons of young, middle-aged, and aged subjects, averaged for each age group (black, blue, and red, respectively, shown as mean ± SEM at each injection level). (C) DRT model output for 4,200 points of the parameter space LHS, shown in 3D projections across the subspace of the excitatory (G EEa and G EEn ) and the inhibitory (G IE ) synaptic weights of pyramidal neurons as their excitability increased (ν ce = 3, 5, 7, and 9 Hz). The numbers in the 3D projection for ν ce = 5 Hz indicate each type of network performance shown in (A), represented as blue, green, yellow, and red dots, respectively, in each graph. (D) Number of points in each LHS maintaining TPA-S as ν ce increased. Dark gray bar indicates the young model cohort. (E) Projections across the (G EEa , G EEn ) subspace showed that points maintaining TPA-S shifted to lower (G EEa , G EEn ) values as ν ce increased. (F) Mean synaptic weights (G EEa , G EEn , G IE , G EIa , G EIn , and G II ) for all points maintaining TPA-S as ν ce increased. G EEa , G EEn , and G IE values shown as purple, pink, and red open circles, respectively; G EIa , G EIn , and G II values shown as dark green, light green, and open blue triangles, respectively. S.E.M. bars lie beneath the symbols.
DRT model (particularly when tuned to the stimulus location), the increased excitability of excitatory neurons was compensated by altering synaptic weights onto those same neurons: with less synaptic excitation and more synaptic inhibition. Yet despite the compensatory synaptic weight changes, fewer networks were able to maintain TPA-S as ν ce increased ( Figure 3D). See also Supplementary Figure 1c for the generalized f -I curve.
These LHS simulations revealed that network function can be maintained somewhat as parameters vary freely to compensate pyramidal neuron hyperexcitability with aging, but the predicted compensation is not consistent with the empirically observed loss of both excitatory and inhibitory synapses on pyramidal neurons with aging. To explore how concomitant changes to firing rates and synapses might affect network function, we first considered all networks maintaining TPA-S when ν ce = 5 Hz as corresponding to a simulated cohort of cognitively healthy young monkeys (the 679 model networks with ν ce = 5 Hz in Figure 3D, called the "young model cohort"). We then examined the performance of the model networks in this cohort on the DRT after perturbing parameters in proportions that matched empirically observed changes with aging (Figure 4). First is the hyperexcitability for pyramidal neurons in vitro described above and in our past work (Chang et al., 2005;Coskren et al., 2015). Second, Peters et al. (2008) reported a ∼10% loss of excitatory and inhibitory synapses in the rhesus dlPFC neuropil by middleage, and a ∼30% loss of both synapse types in aged subjects. This finding is consistent with other studies of synapse and spine loss with aging in the rhesus dlPFC (Hof et al., 2002;Duan et al., 2003;Young et al., 2014).
Accordingly, we induced "middle-aged" and "aged" conditions for pyramidal neuron excitability in the young cohort by increasing ν ce to 7 and 9 Hz, respectively. To model synapse loss, we reduced G EEa , G EEn , and G IE in the young model cohort by 10 and 30%, respectively to create either a middleaged or aged synapse condition. We then induced the two conditions separately and together (forming the "observed data conditions"), and determined whether each transformed model still maintained TPA-S during DRT simulations (Figures 4A-C). A similar synapse loss condition for excitatory synapses was modeled by Wang et al. (2011) to explain an age-related decrease in firing rates during in vivo recordings of the DELAY neurons for monkeys successfully performing the DRT. Our study builds on Wang et al. (2011) by adding inhibitory synapse loss and f -I curve changes with aging, and examining how these alterations affect task performance for a range of parameter values.
DRT performance of the young model cohort was much more sensitive to the excitability aging condition than to the synapse condition (Figures 4A,B). Even under the strongest (aged) synapse condition, over half the networks still maintained TPA-S ( Figure 4A). In contrast, under the excitability condition most of the networks had already lost their capacity for TPA-S by middle age (Figure 4B). In the parameter space exploration above, pyramidal neuron hyperexcitability was partially compensated by decreasing excitatory synapse weights and increasing inhibitory ones. Thus, while Peters et al. (2008) found that inhibitory synapse counts actually decrease with aging, we might still expect a partial recovery of function when the young cohort underwent the excitability and synapse aging conditions simultaneously ("observed data conditions, " Figure 4C). However, this was not the case, as almost none of the middle aged and aged models maintained TPA-S after these perturbations. This does not seem to match reality, since Wang et al. (2011) noted no difficulty in identifying middle-aged and aged monkeys who perform the DRT successfully, and impairment on the more difficult DRSTsp is noticeable but not severe by middle age (Figure 1B).
Since external inputs in vivo are mediated by synapses whose numbers likely also decrease with aging as in Peters et al. (2008), we tested whether reducing the input stimulus current affected the results for the conditions applied in Figures 4A,C  (Supplementary Figures 2a,b). We reduced the stimulus scaling factor, I st , by 10 and 30% for middle aged and aged networks, respectively, in addition to applying the synapse condition (Supplementary Figure 2a) and both observed data conditions (Supplementary Figure 2b). The stimulus reduction did cause a few more of the young cohort models to fail in the middle age and aged conditions, but the overall trend of the results and the number of impaired networks was similar. In particular, DRT performance of the young model cohort was still much more sensitive to the excitability aging condition than to the synapse condition with stimulus reduction (compare Figure 4B and Supplementary Figure 2a).
Recent studies have shown that including short-term synaptic facilitation in discrete and continuous attractor models can reduce the fine-tuning needed to maintain performance (Mongillo et al., 2008;Itskov et al., 2011;Rolls et al., 2013). We added short-term synaptic facilitation to our middle-aged and aged DRT model networks (same parameter values in both groups), as a baseline exploration of how this mechanism might affect DRT performance with aging. The results were qualitatively similar to those without facilitation, except the model cohort networks were affected less by the aging perturbations (Figures 4D-F and Supplementary Figures 2c,d  including stimulus reduction). Under the synapse condition, almost all networks maintained TPA-S both for middle-aged and aged simulations, with even a slight increase from middle-aged to aged when the input stimulus was constant. Networks were simulated time for all sampled networks maintaining TPA-S when ν ce = 5, 7, and 9 Hz (networks shown in Figure 3E). Shown are neurons with preferred directions matching the 0 • stimulus; gray curve shows the mean FR for neurons with a 180 • preferred direction, which fired in none of these networks. Right, mean total synaptic input current to excitatory neurons with the 0 • preferred direction vs. simulation time for the same networks. (B) Mean FR of excitatory neurons for all networks maintaining TPA-S for the model cohort perturbations when the observed data conditions were applied ( Figures 4C,F). Left: without facilitation. Right: with facilitation. (C) Top row: DRT model cohort with applied both observed data conditions and increased excitability of inhibitory neurons. Bottom row: Mean FR of excitatory neurons for all networks maintaining TPA-S. Left: without facilitation. Right: with facilitation. In all panels, young, middle-aged, and aged groups are shown in black, blue, and red, respectively; shaded regions indicate 95% confidence intervals about each mean. still more sensitive to the excitability aging condition than to the synapse condition ( Figure 4D and Supplementary Figure 2c vs. Figure 4E), but under the excitability condition nearly all the middle-aged and about half the aged networks maintained TPA-S. Applying both observed aging conditions together with facilitation led to more networks maintaining TPA-S ( Figure 4F and Supplementary Figure 2d), suggesting that synaptic facilitation can boost DRT performance.
Examining Firing Rates During the Simulated DRT Wang et al. (2011) examined in vivo firing rates of area 46 (dlPFC) pyramidal neurons in rhesus monkeys across the adult age span while they performed the DRT. For monkeys performing the task well, they found the firing rate of DELAY neurons during the cue and delay periods decreased with aging, both for neurons whose preferred and anti-preferred directions matched the cue presented. In our sampling of parameter space (Figure 3) we found that DRT model networks maintaining TPA-S-representing monkeys that perform the DRT successfullyhad lower firing rates during both the cue and delay periods as ν ce increased (Figure 5A, left). This was particularly true for neurons whose preferred direction matched the stimulus, with a substantial decrease from young and middle-aged to aged model networks. This seemed counterintuitive at first, since increasing ν ce makes individual excitatory neurons more excitable in response to a given input current ( Figure 3B). Indeed, hyperexcitability of individual neurons often did lead to over-excitation throughout the network (yellow and red points in Figure 3C). However, the networks that maintained TPA-S as ν ce increased compensated for the increased excitability with a reduction in G EEa and G EEn and increase in G IE (Figure 3F), leading to lower firing rates. Cue period FR in our model was driven much more by the stimulus current than by synaptic weights, and decreased less across the age groups than the delay period FR did. As such, reducing the stimulus current in the sampled middle aged and aged networks as in Supplementary Figure 2 was necessary to match the Wang et al. (2011) cue period results (Materials and Methods). As a whole, these parameter space explorations show a way to unify our in vitro data and in vivo experiments, provided that the total input current to pyramidal neurons during the DRT, I E , is lower in middle-aged and aged monkeys than in young monkeys (Figure 5A, right).
We also examined firing rates in the model cohort simulations perturbed by the observed data conditions with and without facilitation (Figures 4C,F). The main distinction between the parameter space sampling simulations and the cohort simulations was how synaptic inhibition onto the excitatory neurons, G IE , varied as the excitability parameter (ν ce ) increased: G IE increased in the parameter space explorations but decreased in the cohort simulations due to observed inhibitory synapse loss with aging. Indeed, the weakened inhibitory connections in the middle-aged and aged cohorts led to a FR increase during the cue and delay periods of those networks compared to the young cohort (Figure 5B, left), opposing the Wang et al. (2011) results. The increase was even greater when synaptic facilitation was included (Figure 5B, right).
Our young cohort simulations did not match data from the in vitro physiology, synapse counts, and in vivo physiology findings simultaneously, suggesting that one or more age-related changes to the system may be missing from our current model. We simulate one such possibility here, to demonstrate how the current model can help us test hypotheses and incorporate new data as they become available. Specifically, our parameter space explorations showed that an overall increase in inhibition onto the excitatory neurons can compensate for increased firing rates of excitatory neurons. Thus, we examined how network function would change if the observed data conditions were accompanied by a comparable age-related increase of inhibitory neuron firing rates-an empirical question whose answer is not yet known. In models without synaptic facilitation, more middle-aged and aged networks maintained TPA-S when inhibitory neuron firing rates increased with aging ( Figure 4C vs. Figure 5C, top left). Firing rates of excitatory neurons in those simulations also decreased with aging as observed in vivo, consistent with the Wang et al. data (Wang et al., 2011; Figure 5C, bottom left). However, results differed for models that included synaptic facilitation: fewer middle-aged and aged networks maintained TPA-S when inhibitory neuron firing increased (Figure 4F vs. Figure 5C, top right), while (as before) firing rates of excitatory neurons during the DRT increased with aging ( Figure 5C, bottom right). These results suggest that our DRT model might need further refinement to reflect reality better-for example, by also incorporating synaptic depression (Mi et al., 2017) or changes in synaptic plasticity with aging. Regardless, our model allows us to quantify how much the known age-related changes at the single cell level affect cognitive performance, and to predict the extent to which other relevant processes might contribute to working memory impairment with aging.

Pyramidal Neuron Hyperexcitability Led to Reduced Memory Capacity in Simulated DRSTsp Networks
These DRT modeling results add to the wide range of literature using the bump attractor model to study spatial working memory, but the DRT was not administered to the monkeys from this study. Spatial working memory capacity in these monkeys was assessed with the more complex DRSTsp. To see what predictions our model might have for these data, we extended our DRT model to simulate memory retention in a task with three important similarities to the DRSTsp: (1) the network encoded several stimuli simultaneously, with stimuli that are (2) spatially distinct; and (3) presented to the network successively over time. The changes made to extend the DRT model simulations to a task like the DRSTsp are shown in Figure 6A, and a simulated cognitive score (DRST sim ) was assigned to the outcome of each simulation. Examples of model networks with a DRST sim of 4 (perfect performance) down to 1 are shown in Figure 6B. The second network achieved a DRST sim of 3 because it encoded stimuli 1, 2, and 3 correctly but lost persistent activity corresponding to stimulus 3 in the subsequent delay period (red arrowhead). Likewise, the third network (DRST sim = 2) lost persistent activity corresponding to stimulus 2 before stimulus 3 was applied; the final network (DRST sim = 1) encoded stimulus 1 properly but not stimulus 2.
We examined how DRST sim varied across the parameter space of synaptic weights for the young, middle-aged and aged network models (ν ce = 5, 7, and 9 Hz; Figure 6C). Networks achieving a DRST sim of 1 (worst), 2, 3, and 4 (best) are shown, respectively, as blue, yellow, red, and green. As in the DRT, the value of ν ce strongly affected the ability of the DRSTsp model networks to maintain persistent activity tuned to the location of the different stimuli as the number of delay periods and unique stimuli increased. There were many fewer green points in the parameter regime explored, and more blue points, as ν ce increased (As ν ce increased, networks with the best DRST sim performance had lower values of G EEa , G EEn , and G IE than those shown here.). Figure 6D summarizes these counts within our regime, showing a clear leftward shift in the distribution of DRST sim across the parameter space (poorer performance) as ν ce increased. Figure 6E shows how the model networks of Figure 6B were affected by increasing ν ce from 5 (in Figure 6B) to 7 and 9 Hz. In general, as ν ce increased, the initial 0 • stimulus was encoded by a stronger bump (higher firing rate) for the entire DRSTsp simulation, and subsequent stimuli encoded by comparatively weaker bumps. This led to an under-excitation of the network during the second and third delay periods, losing the ability to maintain later stimuli. An example can be seen in the second column of Figures 6B,E. For the young monkey case (Figure 6B, ν ce = 5 Hz), the maximum firing rate of the first bump was around 75 Hz; this and the width of the first bump decreased slightly as the simulation continued. Also, as noted above, persistent activity for stimulus 3 was lost during the third delay period. For the middle-aged case (Figure 6E top row, ν ce = 7 Hz), the maximum firing rate of the first bump maintained a maximum firing rate near 150 Hz and a similar width throughout the simulation, but the second bump was lost early in the subsequent delay period. This trend was even stronger for the aged case (Figure 6E bottom row, ν ce = 9 Hz), where the first bump remained wide with a maximum firing rate of 250 Hz throughout the simulation but none of the subsequent stimuli were encoded.
To see how synaptic parameters compensated to maintain DRSTsp performance as ν ce increased, we looked at the mean values of the synaptic weights for all networks with perfect task performance (DRST sim = 4). Values of both the excitatory and inhibitory synaptic weights for excitatory neurons were lower as ν ce increased ( Figure 6F top; standard error bars lie beneath the symbols). Excitatory synaptic weights for inhibitory neurons decreased and inhibitory synaptic weights increased as FIGURE 6 | Memory retention model of the DRSTsp shows impairment under neuronal hyperexcitability. The ability to maintain the memory of a successively increasing number of stimuli strongly depended on the excitability of individual pyramidal neurons, but changes in the synaptic weighs partially compensated for the increased neuronal excitability. (A) Simulated task included a pre-cue period, with four cue periods separated by three delay periods. Green circles indicate the stimuli presented in each cue period. (B) Examples of DRSTsp model output, with FR shown in color (log scale) for neurons labeled by their preferred direction vs. simulated time. From left to right: DRST sim = 4, optimal performance, all four stimuli were encoded correctly and recalled; DRST sim = 3, where the first three stimuli were encoded correctly, but persistent activity representing the third was lost during the subsequent delay period; and DRST sim = 2 and DRST sim = 1 where only the first two or one stimuli were encoded correctly and recalled. Cue periods shown between white lines; arrowheads indicate when persistent activity was lost. (C) DRSTsp model output for 4,200 points of the parameter space LHS, shown in 3D projections across the (G EEa , G EEn , G IE ) subspace as ν ce increased. Points with DRST sim = 4, 3, 2, and 1 are shown as green, red, yellow, and blue dots, respectively. (D) Histogram of DRST sim values in the LHS as ν ce increased (ν ce = 5, 7, and 9 Hz shown as black, blue, and red bars, respectively). (E) The first three simulations from (B), repeated for ν ce = 7 (top row) and 9 Hz (bottom row). (F) Mean synaptic weights for all LHS points with DRST sim = 4 for ν ce = 5, 7, and 9 Hz. Top: G EEa , G EEn , and G IE shown as purple, pink, and red open circles. Bottom: G EIa , G EIn , and G II shown as dark green, light green, and open blue triangles. S.E.M. bars lie beneath the symbols. ν ce increased ( Figure 6F bottom; standard error bars lie beneath the symbols). Overall, for the DRSTsp model, both panels showed that synaptic excitation as well as inhibition to excitatory neurons decreased to partially compensate for the increased excitability of individual excitatory neurons.
"Aging" the DRSTsp Model Networks Led to Working Memory Impairment as Observed Empirically As with the DRT model, we defined a simulated "young monkey cohort" for the DRSTsp model and "aged" the model networks in this cohort by applying increased excitability and synaptic loss conditions (Figure 7), choosing 40 levels of parameter perturbations throughout the simulated age span. We induced the two observed data conditions separately and together-and to match Figure 5C top right, also increased the excitability of inhibitory neurons (Since reducing the stimulus current had a minimal effect on performance during the synapse condition perturbations of the DRT, we kept the stimulus current constant across all DRSTsp simulations.). We then computed the corresponding DRST sim for each transformed member of the young cohort. To summarize the change in performance across the simulated age span, we calculated the relative change in the mean DRST sim from the young to the aged models under each perturbation condition. This facilitated comparison between performance on the simulated task (four stimuli spaced evenly around a ring), and the real DRSTsp (18 possible stimuli, organized spatially into three rows of six wells each).
All four cases showed a strong linear relationship between DRST sim and parameter perturbations. As in the DRT model above, the DRSTsp performance of the young model cohort was much more sensitive to the excitability condition ( Figure 7B) than to the synapse condition ( Figure 7A). While the synapse condition showed a 4% improvement with aging for the DRSTsp model (R 2 = 0.98, p < 0.0001), the excitabilityonly condition showed a 45% decrease in DRST sim with aging (R 2 = 0.98, p < 0.0001). Applying the two observed conditions together led to a 14% impairment with aging ( Figure 7C; R 2 = 0.82, p < 0.0001); applying the observed data conditions plus increasing the excitability of inhibitory neurons led to a 17% impairment ( Figure 7D; R 2 = 0.81, p < 0.0001). Thus, both cases of applying the observed data conditions were comparable to the 20% impairment exhibited by the empirical data ( Figure 1B). These results are also consistent with Figure 6F: the empirically-observed decrease in synaptic excitation and inhibition to excitatory neurons partially compensated the increased excitability of the excitatory neurons, so that DRST sim impairment with aging was less severe. The model also suggests that hyperexcitability of inhibitory neurons would impact performance on the DRT more than it would the DRSTsp.

DISCUSSION
This study used computational modeling to explore how single pyramidal neuron hyperexcitability and excitatory and inhibitory synapse loss observed with aging in layer 3 of the rhesus monkey dlPFC might contribute to spatial working memory impairment. First, we presented new empirical data from rhesus monkeys across the adult age span, consistent with past studies, showing that working memory impairment and increased FRs of neurons occur in some middle-aged monkeys. This replicated previous studies showing a decline in cognitive performance from young to aged monkeys (Herndon et al., 1997;Moore et al., 2003Moore et al., , 2006Konar et al., 2016;Motley et al., 2018), and a concomitant increase in R n and FR of dlPFC pyramidal neurons (Chang et al., 2005;Coskren et al., 2015). These data motivated the modeling, and were used to constrain a bump attractor network model of the DRT (Wimmer et al., 2014) and our network model of working memory retention in the DRSTsp. These models predict that the observed age-related firing rate increase and synapse loss in rhesus dlPFC (Peters et al., 2008) are sufficient to induce significant impairment of spatial working memory function.
There is evidence that persistent neural activity during delay periods of a working memory task is important for encoding the memory of a stimulus (reviewed in . When age-related alterations to the excitability properties of neurons comprising the working memory network occur, the network's capacity for maintaining persistent activity-keeping activity tuned to the original stimulus-is impaired (Wang et al., 2011). To provide insight into the relationship between neuronal firing rate and synaptic changes seen in aging and working memory decline, we examined the parameter space of a bump attractor model of the DRT, then extended the model to simulate how aging might affect memory retention of multiple stimuli in a more complex task-the DRSTsp. The region of parameter space that simulated the DRT successfully was largest when the pyramidal neuron excitability parameter (ν ce ) was tuned to fit the young FR data. The network model had its own built-in compensatory mechanisms: when pyramidal neuron excitability increased, the region with successful DRT models shifted to lower excitatory and higher inhibitory synapse weights. This tradeoff allowed networks to control the overall excitatory:inhibitory balance needed to maintain TPA. These compensations were consistent with empirical data showing an age-related decrease in excitatory synapses, but not consistent with the reported decrease in inhibitory synapses (Peters et al., 2008). However, this compensatory mechanism led to lower firing rates during the delay period, as well as during the cue period when the stimulus current was also reduced, unifying our in vitro observations with the in vivo data (Wang et al., 2011). We predict that the network compensates higher FR in individual pyramidal neurons with aging by lower total synaptic input current to each neuron, and that the stimulus current could also be decreasing with aging.
To test how hyperexcitability of pyramidal neurons and loss of excitatory and inhibitory synapses with aging affected network output, we ran a second set of simulations: "aging" model neurons by decreasing both excitatory and inhibitory synaptic strengths and increasing the FR of the excitatory neurons. DRT model performance (represented by the number of networks performing the task successfully) was substantially impaired in middle-aged and aged simulations, more than reported in a sample of rhesus monkeys (Wang et al., 2011), and the FR increase affected DRT performance much more than the loss of synapses did. Adding synaptic facilitation restored successful function to many middle aged and aged DRT models. However, only one of the aging cohort scenarios we tested captured the hyperexcitable pyramidal neuron f -I curves observed in vitro, dlPFC synapse loss observed with electron microscopy, and the in vivo FR decrease during the DRT cue and delay periods from Wang et al. (2011): the one without synaptic facilitation but added hyperexcitability of inhibitory neurons to restore excitatory:inhibitory balance. One next step is to examine the excitability of dlPFC interneurons with aging, but our simulations raise several other questions about the DRT. Is there an age-related reduction in overall input to dlPFC during the cue period? How central a role does synaptic facilitation play during the task? How do age-related changes in synaptic plasticity, discussed below, affect DRT performance? Might neuromodulation (Arnsten et al., 2012;Davis et al., 2017;Wang et al., 2019) during working memory task performance reduce the pyramidal neuron hyperexcitability observed in vitro? Computational models provide an essential means for testing and refining hypotheses as future data become available.
Finally, we extended the DRT bump attractor model to create a model of working memory retention in a task like the DRSTsp. Other network models of working memory have simulated several stimuli at once (Edin et al., 2009;Wei et al., 2012;Mi et al., 2017). Yet, while the DRSTsp has been used for many years to evaluate spatial working memory capacity in rhesus monkeys (Herndon et al., 1997;Luebke et al., 2004;Chang et al., 2005;Moore et al., 2005Moore et al., , 2017Peters et al., 2008;Luebke and Amatrudo, 2012)-including in the present study-to our knowledge this is the first computational model of this task. Increasing pyramidal neuron FR strengthened the response to the initial spatial cue and reduced the ability of the DRSTsp model networks to encode memories of subsequent spatial cues, so that middle-aged and aged networks had lower DRST sim than young simulations. As with the DRT model, "aging" the young DRSTsp model cohort led to severe impairment when increasing pyramidal neuron excitability alone, but not synapse loss alone. The maintenance of the DRSTsp model performance under synapse loss alone is consistent with Peters et al. (2008), which found no correlation between the numerical density of synapses and DRSTsp impairment. Simultaneously varying synapse loss and pyramidal neuron hyperexcitability within ranges observed with aging, both with and without hyperexcitability of interneurons, the amount of DRST sim impairment was similar to that seen in DRSTsp scores across the adult rhesus monkey life span. Thus, our model predicts that pyramidal neuron hyperexcitability and synapse loss may be sufficient to explain empirically observed levels of spatial working memory impairment. Future experiments will examine whether the amount of synapse loss reported in Peters et al. (2008), is also seen in the dlPFC of a subset of monkeys from this study, and whether there are correlations between those data, AP firing rates, and spontaneous postsynaptic currents.
There is an open debate about the role of synaptic mechanisms for working memory Lundqvist et al., 2018). While delay activity is recognized as important to working memory, a recent perspective posited that sparse spiking activity and synaptic plasticity between spike times-rather than asynchronous persistent activity-might be the mechanism actually responsible for memory maintenance . This view is supported by several modeling studies (Sandberg et al., 2003;Mongillo et al., 2008;Lundqvist et al., 2011). Our models lend further support to an essential role for both persistent activity and synaptic facilitation for working memory tasks involving multiple cues. Synaptic plasticity (both facilitation and depression) help overcome limitations of continuous attractor models, such as filtering out distractors (Compte et al., 2000) and remembering multiple items (Itskov et al., 2011;Wei et al., 2012;Hansel and Mato, 2013;Rolls et al., 2013;Mi et al., 2017) as in our DRSTsp model. Bump attractor models of visuo-spatial working memory have shown that short-term synaptic facilitation in recurrent connections of excitatory neurons slows down the drift of the bump often seen in continuous attractor models (Itskov et al., 2011). Hansel and Mato (2013) showed that short-term synaptic facilitation might lead to the direction selectivity of synaptic weights and network bistability in a DRT model. In vitro experiments have shown synaptic facilitation in PFC pyramidal neurons (Hempel et al., 2000;Wang et al., 2006). Several other studies have shown age-related changes in synaptic plasticity (not modeled here), some correlating these changes with learning and memory deficits (Norris et al., 1998a,b;Bach et al., 1999;Rosenzweig and Barnes, 2003;Thibault et al., 2007Thibault et al., , 2013Gant et al., 2011Gant et al., , 2014Gant et al., , 2015. This network-level study represents an extension of our past modeling studies at the single-neuron level on parameters affecting firing properties of individual pyramidal neurons with aging (Coskren et al., 2015;Luebke et al., 2015;Rumbell et al., 2016), in different cortical regions , and in neurodegeneration (Goodliffe et al., 2018). The elegance of these network models lies in their ability to represent empirically observed changes in FR and synaptic weights with aging as perturbations to a small number of network parameters. Based on our findings with two f -I curves, replacing our firing rate model neurons with a spiking neuron model (Stein, 1967;Hodgkin and Huxley, 1990;Izhikevich, 2003;Teeter et al., 2018) in our DRT and DRSTsp networks should yield qualitatively similar results. Our future network modeling will incorporate age-related changes in passive and active channel parameters predicted in our and other past studies (Wang et al., 2011;Coskren et al., 2015;Rumbell et al., 2016), and other empirically observed changes in aging dlPFC including increased slow after hyperpolarization current ; white matter changes (review: Luebke et al., 2010;Kubicki et al., 2019); and dysregulation of calcium homeostasis (Foster, 2007;Toescu and Vreugdenhil, 2010;Oliveira and Bading, 2011). Age-related changes to synaptic strength and neurotransmitters reduced the probability of short-term memory recall from a discrete attractor network model (Rolls and Deco, 2015); similar perturbations in a spiking-network version of our DRT and DRSTsp models would be extremely insightful.
The main limitation of our current DRSTsp model-assuming that the network correctly identifies the novel stimulus among several retained-is also a limitation of the classical bump attractor as a whole. Experiments have shown that decisionmaking and attentional tasks involve connections between dlPFC and other cortical regions including posterior parietal cortex, anterior cingulate cortex, and secondary motor areas (Medalla and Barbas, 2010;Oemisch et al., 2015;Marvel et al., 2019;Tsunada et al., 2019). These connections may represent different aspects of a task like DRSTsp, such as the monkey's internal model of itself and the task; its level of attention; signals of reward and error; and motor traces that might correspond to rehearsal strategies between active stimulus stages of the DRSTsp. Recent theoretical studies have shown that model neurons with mixed selectivity, responding to sensory stimuli as well as internal states of the system, can perform complex cognitive tasks (Maass et al., 2002;Salinas, 2004;Rigotti et al., 2010;Chaisangmongkon et al., 2017). Such units could be well-suited for extending the DRSTsp model introduced here to simulate choice of the novel spatial cue and related internal signals, plus the full 3 × 6 array of wells presented to the monkeys during the DRSTsp. There is evidence that aging affects information flow between cortical regions, which can contribute to working memory impairment (Engle et al., 2016;Lee et al., 2016;Proskovec et al., 2016;King et al., 2018;Koen et al., 2019). Any such age effects must be examined in concert with the changes to dlPFC synapses and neuronal excitability, to predict the extent to which these concomitant changes compound vs. partially compensate each other.
In summary, our computational models of two spatial memory tasks predicted that empirically observed changes in pyramidal neuron excitability and synapse loss lead to spatial working memory impairment with aging. Computational studies such as this provide a means to evaluate connections between physiological and behavioral-and in vitro vs. in vivo-experiments, extending the impact provided by animal studies.

ETHICS STATEMENT
The animal study was reviewed and approved by Boston University Institutional Animal Care and Use Committee (IACUC).

AUTHOR CONTRIBUTIONS
SI developed and implemented the models, analyzed the results, carried out the statistical analysis of the data, and wrote the manuscript. JL performed the experiments, wrote the manuscript, and designed the study. WC performed the experiments. DD designed the statistical analysis of the data. CW carried out the statistical analysis of the data, wrote the paper, and designed the study.