Neurodegeneration exposes firing rate dependent effects on oscillation dynamics in computational neural networks

Traumatic brain injury (TBI) can lead to neurodegeneration in the injured circuitry, either through primary structural damage to the neuron or secondary effects that disrupt key cellular processes. Moreover, traumatic injuries can preferentially impact subpopulations of neurons, but the functional network effects of these targeted degeneration profiles remain unclear. Although isolating the consequences of complex injury dynamics and long-term recovery of the circuit can be difficult to control experimentally, computational networks can be a powerful tool to analyze the consequences of injury. Here, we use the Izhikevich spiking neuron model to create networks representative of cortical tissue. After an initial settling period with spike-timing-dependent plasticity (STDP), networks developed rhythmic oscillations similar to those seen in vivo. As neurons were sequentially removed from the network, population activity rate and oscillation dynamics were significantly reduced. In a successive period of network restructuring with STDP, network activity levels were returned to baseline for some injury levels and oscillation dynamics significantly improved. We next explored the role that specific neurons have in the creation and termination of oscillation dynamics. We determined that oscillations initiate from activation of low firing rate neurons with limited structural inputs. To terminate oscillations, high activity excitatory neurons with strong input connectivity activate downstream inhibitory circuitry. Finally, we confirm the excitatory neuron population role through targeted neurodegeneration. These results suggest targeted neurodegeneration can play a key role in the oscillation dynamics after injury. Author Summary In this study, we study the impact of neuronal degeneration – a process that commonly occurs after traumatic injury and neurodegenerative disease – on the neuronal dynamics in a cortical network. We create computational models of neural networks and include spike timing plasticity to alter the synaptic strength among connections as networks remodel after simulated injury. We find that spike-timing dependent plasticity helps recover the neural dynamics of an injured microcircuit, but it frequently cannot recover the original oscillation dynamics in an uninjured network. In addition, we find that selectively injuring excitatory neurons with the highest firing rate reduced the neuronal oscillations in a circuit much more than either random deletion or the removing neurons with the lowest firing rate. In all, these data suggest (a) plasticity reduces the consequences of neurodegeneration and (b) losing the most active neurons in the network has the most adverse effect on neural oscillations.

injury dynamics and long-term recovery of the circuit can be difficult to control experimentally, 23 computational networks can be a powerful tool to analyze the consequences of injury. Here, we use the 24 Izhikevich spiking neuron model to create networks representative of cortical tissue. After an initial 25 settling period with spike-timing-dependent plasticity (STDP), networks developed rhythmic oscillations 26 similar to those seen in vivo. As neurons were sequentially removed from the network, population activity 27 rate and oscillation dynamics were significantly reduced. In a successive period of network restructuring 28 with STDP, network activity levels were returned to baseline for some injury levels and oscillation 29 dynamics significantly improved. We next explored the role that specific neurons have in the creation and 30 termination of oscillation dynamics. We determined that oscillations initiate from activation of low firing 31 rate neurons with limited structural inputs. To terminate oscillations, high activity excitatory neurons with 32 strong input connectivity activate downstream inhibitory circuitry. Finally, we confirm the excitatory 33 neuron population role through targeted neurodegeneration. These results suggest targeted 34 neurodegeneration can play a key role in the oscillation dynamics after injury. 35 36 Introduction 49 Traumatic Brain Injury (TBI) is a prominent cause of disability in the US (1). Perhaps because of growing 50 awareness of the consequences of TBI, emergency department visits for TBI have increased 47% from 51 2007 to 2013 (2). Although many of these injuries produce no long-term deficits, a fraction of injuries 52 produce cognitive and psychological impairments that can last years after the original insult (3-5). As a 53 result, more than 5 million people in the US live with significant consequences of TBI, contributing to an 54 estimated 70 billion dollars annually in medical and non-medical costs (6). Effective recovery from TBI 55 remains a challenge because no two injuries are exactly alike, leading to a unique injury and recovery 56 pattern for each TBI patient. 57 One key feature in TBI recovery is how the structural and functional networks in the brain evolve over 58 time after injury to guide the cognitive recovery processes (7)(8)(9)(10)(11). Recent studies have shown 59 alterations in brain circuitry after moderate and severe injury affect the coordination among functional 60 brain networks (12). With the development of models to predict the overall changes in brain networks 61 during different tasks, there is an emerging consensus that the dynamic network that connects different 62 brain regions can influence cognitive and psychological alterations after injury (13-16). However, the 63 unique circumstances that cause each TBI make it difficult to predict which injuries will likely lead to 64 long-term changes in brain function. These challenges exist especially at the cellular scale, where the 65 neuronal degeneration that can occur days to weeks after a TBI can affect the function of local 66 microcircuits throughout the brain (17,18). 67 To this end, computational models can be a useful tool to understanding how neuronal damage can 68 ultimately contribute to the impairments in circuit function after TBI. In general, these models can 69 account for the mechanisms of acute injury to the network (e.g., primary axotomy, membrane 70 permeability changes, receptor dysfunction) and secondary changes that can also trigger neuronal loss 71 These changes were not sufficient to return to baseline, but partially restored function. Dotted vertical 131 line at 70% and 75% injury denote the last injury that all simulations contained oscillations in pre-and 132 post-plasticity simulations respectively. (E) FWHM differences after plasticity were significantly different 133 from baseline and significantly different from immediately post-injury in the range of 40-70% injury. 134 135 After observing that spike-timing dependent plasticity helped introduce resilience to damage, we next 136 explored if there were specific connectivity features of individual neurons that influenced or explained 137 part of this resilience. For each neuron in an undamaged network, we computed a neuron connectivity 138 index as the normalized difference of total synaptic input strength and the total output strength. This 139 neuron-connectivity index correlated with the average neuronal firing rate; neurons with high index 140 showed higher firing rates than neurons with low index ( Figure 3A). The relationship between 141 input/output strength and firing rate was stable while allowing the synaptic weights to adjust via STDP 142 over 4 simulation hours. In addition, these relationships did not change across a broad range of neuronal 143 network parameters that included synaptic strength, neuron parameters (a-d) that could change 144 neuronal type (33), and connection number among neurons in the network. 145 We then sought to find if either the initiation or termination of the oscillations were related to neuronal 146 activity and, in turn, connectivity strength. Using our definition of the beginning and end of an oscillation 147 (see Methods; representative oscillation appears in Figure 3B, C), we divided an oscillation period into 148 quintiles. In general, excitatory neurons with low inputs relative to their outputs (i.e., low index) fired 149 primarily during the initiation period of an oscillation, while high input strength excitatory neurons fired 150 near the peak of an oscillation and activated inhibitory neurons to arrest the oscillation. The average 151 firing rate of excitatory neurons within each of the first four quintiles significantly increased, while the 152 firing rate significantly decreased in the fifth quintile (One way ANOVA with Tukey Kramer post-hoc 153 (p<.001; Figure 3D)). In comparison, we could not identify any dependence on firing time and oscillation 154 period for inhibitory neurons ( Figure 3D). Excitatory neurons most commonly fired during the peak of 155 the oscillation, decreasing on both sides of peak ( Figure 3E). Inhibitory neurons had a preference to fire 156 later in the oscillation, peaking between the 3 rd and 4 th quintile ( Figure 3E). Together, these results 157 indicate that the initiation of an oscillation corresponds with the activations of neurons with low firing 158 rates, while the termination of the oscillation begins with the simultaneous recruitment of neurons with 159 high firing rates and activating a significant fraction of the inhibitory neurons. 160 peak oscillation magnitude, with decreasing activation after the peak. Inhibitory neurons had delayed 173 activation, primarily firing at peak or late in the oscillation. 174 175 As excitatory neurons with specific connectivity and activity patterns correlated with oscillation 176 dynamics, we next considered if the targeted deletion of either activity type would alter the neural 177 circuit dynamics. Using the random deletion of neurons as a comparison, we explored the change in 178 functional network characteristics that would occur if we progressively deleted the neurons with the 179 lowest average activity. With this strategy, we found that average activity rate would significantly 180 decrease when damage exceeded 5% of the network ( Figure 4A-B). Similar to random deletion, less 181 damage was required to significantly changes oscillation frequency ( Figure 4C; damage > 5%) than 182 activity rate. The deletion of the lowest firing rate neurons additionally led to significant change in the 183 average width of an oscillation ( Figure 4E). Plasticity returned the average activity in the network to 184 baseline up to 15% damage; beyond this level, plasticity did improve average firing rate up to 90% injury 185 ( Figure 4A). In comparison, plasticity improved the oscillation frequency (and magnitude) over a range of 186 the injury levels (10-65% and 5-55% respectively) but did not reach baseline levels. 187

Figure 4. Effect of deleting relatively inactive neurons is partially recovered with plasticity (A) 189
Representative raster plots before injury, at 50% low firing rate (LFR) excitatory neuron injury, and after 190 network restructuring with STDP . (B) LFR damage produced a rapid decline to baseline firing after injury 191 that significantly recovered with STDP. (C,D) This restoration was primarily due to a recovery of network 192 oscillation frequency and magnitude. (E) Changes in oscillation width occurred relative to baseline, 193 primarily following damage exceeding the threshold of full firing rate restoration. 194 195 In contrast to these results, the progressive deletion of neurons with the highest activity rate did not 196 significantly change the overall average activity of the network until more than 20% of neurons were 197 removed ( Figure 5A-B). Removing the neurons with highest activity led to a progressive decrease in 198 oscillation frequency that was significantly different than deleting the same fraction of low firing rate 199 neurons (t-test with Bonferonni correction for multiple comparisons; p<.001; Figure  Here, we modeled the effects of neurodegeneration on the functional dynamics of neural circuitry. We 214 utilized an established spiking neural network model to investigate how damage and plasticity impact 215 network firing and oscillatory behavior. We show that neurodegeneration, regardless of whether it was 216 random or based on neuronal firing rate, significantly decreased network activity and oscillation 217 dynamics. For all types of neurodegeneration, activity was significantly restored with spike-timing 218 dependent plasticity. Deleting highly active neurons led to a marked increase in oscillation duration, 219 while deleting neurons with low activity reduced the initiation of oscillations without decreasing their 220 frequency. These results suggest that the degeneration or inactivation of specific neuron activity profiles 221 can differentially affect oscillation dynamics of neuronal circuits. 222 We used several simplifying assumptions to examine the potential impact of neurodegeneration. First, 223 we use generalized topologies and neuron spiking behaviors of cortical neurons based on simplistic rules 224 developed from observations in vivo (28,32,33,35,36). We recognize that neurons within specific brain 225 regions can vary greatly in the connectivity preferences and spiking behavior (37-39), and therefore 226 some of our observations on neuronal degeneration may not apply universally to all brain regions. We 227 did explore if our observations were influenced by different distance-dependent connection algorithms 228 and found that the general decline in activity and the restorative effect of plasticity did not depend on 229 the initial spatial connections in the network. Although we anticipate that our general findings could 230 inform predictions for networks of more complex, region specific topologies, many of the regions 231 commonly damaged in TBI (e.g., hippocampus, thalamus, and cingulate (29,30)) do not have clear 232 estimates of neuronal connectivity. Once available, we envision creating specific computational network 233 models to assess the dynamics would differ with deletion in these specific regions. 234 Second, we used Izhikevich integrate-and-fire neuron model to approximate neuron activity, which may 235 not fully represent the complexity of real circuitry. Because the Izhikevich model accurately predicts 236 neuronal spike times for a variety of neuron spiking behaviors (33,40,41), this model should be 237 representative of in vivo neural activity for the metrics we utilize in this study. Extending the current 238 work into a more computationally complex neuron model with higher-order biological features 239 (reviewed in (42)) could provide additional information into the timing of activation and synaptic inputs 240 across an individual neurons (43). However, these changes would likely affect all neurons in the network 241 similarly and will not significantly impact the overall estimates of activity and oscillatory behavior. 242 Finally, we considered whether the broad distribution of neuronal types employed in larger scale 243 Izikevich integrate-and-fire models were important to include in this study. However, without any 244 insight into how more homogeneous circuits remodel in response to trauma, we believed this additional 245 complexity would prevent one from gaining clear insight into which types of neuronal deletion would be 246 particularly damaging to circuit dynamics. 247 In general, our observed neural activity patterns match the general spiking patterns and network wide 248 oscillations in computational networks of similar size (e.g., (27,33)), and even more complex networks of 249 much larger size and complexity (32). One key characteristic we observed in our simulations was the 250 coordinated activation of a subset of neurons in regular periodic intervals over the entire simulation 251 period. Again, these waves of neuronal activation are distinct from a near simultaneous activation or 252 bursting of the network that can appear in some studies (27). These oscillations of neuronal activity, 253 where 10-15% of the network was activated, were the first to show a significant decrease in oscillation 254 frequency after neurons were deleted from the network. Although neuron characterization traditionally 255 depends on functional spiking properties of neurons (44-46), our results show that plasticity is an 256 important mechanism to produce functional diversity. Moreover, our results identified subpopulations 257 of neurons that could either help trigger or suppress these periods of high network activity. To our 258 knowledge, showing that neurons with low activity rates preferentially activate oscillations in a network 259 has not been reported in past modeling studies, nor are we aware of past reports showing that neurons 260 with high firing rates are important for quieting periods of high activity. Together, these data suggested 261 that targeting of specific neuronal populations for degeneration would preferentially affect neuronal 262 dynamics, a prediction we confirmed with subsequent simulations. 263 The rhythmic oscillations may also have important consequences on the synchronization, or coherence, 264 of activity across brain regions. Coherence among neuron populations has been shown to be important for attention and memory (47-49), cognitive processes that are commonly affected after traumatic 266 brain injury. Selective degeneration of low firing rate neurons reduces the likelihood of initiating an 267 oscillation, in turn lowering the coherence with other brain regions downstream of the injured 268 microcircuit. Therefore, losing this neuronal subpopulation would appear to play a significant role 269 affecting information relayed across brain regions, an aspect of information processing that has 270 appeared in network-based studies of TBI (13,50). In comparison, losing neurons with high firing rates 271 would appear to have less consequence on the coordinated oscillations among regions because these 272 neurons would not affect the emergence of an oscillation, and only slightly lengthen the duration of an 273 oscillation. However, we cannot completely discount the impact of losing high firing rate neurons 274 because the lengthening of a specific oscillation may impede the propagation of sequential information 275 across nodes in a network. Together, these point to a potential role for a small number of neurons to 276 play a larger role in relaying information, via oscillations, across several interconnected microcircuits. 277 Our general finding that spike-timing dependent plasticity (STDP) is a key mechanism to re-stabilize 278 network dynamics provides a potentially new role for STDP in the injured brain. As a primary mechanism 279 associated with Hebbian learning, STDP is typically considered as a mechanism to restructure the 280 synaptic connections in a network after a training stimulus (51,52). The return of activity to a damaged 281 network using STDP is reminiscent of how homeostatic plasticity allows a healthy network to gravitate 282 towards a target activity rate (53). Interestingly, homeostatic plasticity may play a very different and 283 destabilizing role in networks after either focal or more diffuse deafferentation of neurons (27), leading 284 these networks into brief bursts of activity the resemble interictal discharges that appear in 285 posttraumatic epilepsy (54,55). However, this rebalancing of networks with STDP has its limits, as the 286 relative success of recovering initial dynamics is not complete at the highest injury levels. In light of 287 several reports showing that one form of plasticity -long term potentiation (LTP) -is lost after 288 traumatic injury in vivo and in vitro (56,57), our results emphasize the importance of therapeutic 289 strategies to help promote plasticity after injury and regain initial network dynamics. Similarly, given 290 that oscillations can be important to establish coherence among brain regions, steps to maintain 291 plasticity in an injured network would likely improve network communications in the traumatically 292 injured brain. 293 Our results also provide some suggestions on how local damage in the brain may affect both local and 294 global brain dynamics. Cognitive disruptions from TBI are frequently viewed as changes in the network 295 structure among regions in the brain (reviewed in (13,50,58-60). Commonly, past studies focus on the 296 functional and structural deficits in the connections among nodes in the network resulting from diffuse 297 axonal injury (29,50,61). With new techniques in medical imaging and network theory, we can begin to 298 understand how changes in the connectivity between regions can impact higher level cognitive function 299 (62-64). However, these approaches rely on maintaining function at the node (microcircuit) level after 300 injury. Certainly, we know there are regions of the brain that can impart large cognitive deficits simply 301 with their own malfunction (50,65). From the current work, we know that neurodegeneration can 302 impact both network activity and neural oscillations in the node. In combination, our results suggest 303 that damage to one node (microcircuit) could indirectly influence the coordination of activity across 304 many connected brain areas. 305 It is certainly plausible to explore some of the unique consequences of local neurodegeneration in 306 broader brain networks using oscillator or neural mass models to link structural and functional networks 307 of the brain (66,67). Although these models capture gross behavior of network dynamics, the current 308 formulation of these models lack the nodal accuracy to determine how perturbations caused by injury 309 can impact the larger network function. Similar to work showing how local gamma activity could create 310 biologically realistic BOLD correlations (32), changes in network oscillation frequency or oscillation width 311 would likely have broader impact beyond the local network. Particularly important would be examining the implications of an inconsistent oscillatory rate, as seen in our model, in a connected oscillatory 313 model. To our knowledge, this feature is not commonly explored in neural mass or oscillator-based 314 models. However, there is evidence to suggest that these models have synchronization properties that 315 are sensitive to nodal dynamic changes, but oscillators with individually variable frequencies have yet to 316 be investigated (68,69). Similarly, both neural mass and oscillator-based models can potentially benefit 317 from further analysis of biological plasticity mechanisms to repair from damage, given the role we found 318 plasticity in stabilizing or practically recovering nodal dynamics. 319 Overall, this study indicates that neurodegeneration alters population-level activity and network 320 oscillations, with subpopulation-dependent changes to oscillation frequency or duration. These changes 321 in network dynamics can be significantly recovered with spike-timing dependent plasticity. We 322 anticipate that future work brain network dynamics to specific patterns of damage will be valuable for 323 gaining insight into specific patterns of damage causing longlasting alterations in brain dynamics, and 324 other patterns of damage that may produce only temporary changes in neural dynamics. At a higher 325 level, discriminating between these two injury patterns can help identify injuries that could cause lasting 326 cognitive deficits much earlier than currently possible, pointing to an opportunity to treat and improve 327 outcome in a vulnerable population of TBI survivors. 328

Modeling a representative cortical circuit 330
To investigate the connection between neural degeneration and functional network activity, we 331 constructed computational neural networks of integrate-and-fire neurons ((33); Summary in Figure 6). 332 Networks of 1000 neurons were constructed with 80% regular-spiking, excitatory neurons and 20% fast-333 spiking inhibitory to mimic the ratios commonly used to model cortical circuits (28,32,33). These 334 neurons follow a system of ordinary differential equations to track membrane potential, membrane 335 recovery, and threshold-based spiking as follows: 336 ′ = .04 2 + 5 + 140 − + 337 Where v represents the membrane potential in millivolts, and u is the membrane recovery variable. 340 Parameters a, b, c, and d were set to create heterogeneous regular-spiking excitatory neurons and fast-341 spiking, low-threshold inhibitory neurons as in Izhikevich 2003 ( Figure 6B). We allowed parameters a, b, 342 c, and d to vary within a tight range for each neuron to avoid network behavior resulting from a 343 homogeneous neuron population. We used a fixed timestep (0.2 millisecond) and a forward Euler 344 method to compute v and u over time. 345 Consistent with average firing rates in vivo and previous models using Izhikevich neurons (28,35,36), we 346 used a gamma distribution (k, θ = 2, ½) to randomly (f= 1 Hz) inject currents into individual neurons 347 within the network. This stimulation was strong enough to cause the neuron to fire and send AMPA-or 348 GABA-based synaptic signals to downstream targets. Synaptic currents were modeled as exponential 349 decays from AMPA or GABA receptors, with τ = 5 msec. AMPA and GABA receptors were initially set to 350 create EPSP/IPSP in line with past in vivo recordings (70). Repeated input stimuli to were attenuated at 351 40% immediately after a spike occurred (τ=150 msec) to model desensitization in the neuron 352 population. 353 To avoid bias from edge effects in seeding neuron position, we placed neurons randomly on the surface 354 of a unit sphere ((28); Figure 6C). The number of outputs for each neuron, which was drawn from a 355 normal distribution with an average of 100 total outputs and inputs per neuron, varied slightly from each neuron to mimic features estimated in cortical circuits (10% variance; (71)). Neurons were 357 randomly connected to each other across the surface, producing network properties of a classic  Renyi random graph (72). In a subset of simulations, we examined the effects of weak and strong 359 distance-dependent connections and found that our main findings were unchanged. As a result, we 360 present only the simulations using a random connection topology. Finally, we implemented synaptic 361 transmission delays that were proportional to the distance between two neurons along the arclength of 362 the sphere and set 8 ms as the maximum delay, consistent with in vivo recordings (73-77). 363 In neural networks, synaptic connection strengths adapt according to different models of synaptic 364 plasticity. Among the models available, we chose to implement spike-timing dependent plasticity (STDP) 365 because of its critical role in learning and the potential role this feature may play in cognitive deficits 366 after traumatic injury (56,78,79). We used the Song model of STDP (34,80), in which the synaptic 367 strength was adjusted based on the relative timing of synaptic inputs to a neuron and the subsequent 368 action potential firing of the target neuron. Mathematically, this can be described as: 369 Where w is the synaptic weight of the connection between the pre-and postsynaptic neuron, A+ and A-,

374
Given that this formulation of STDP leads to a bimodal distribution of synaptic weights, in which weights 375 approach either the minimum or maximum possible weight (80), we seeded the synaptic strengths in 376 our networks using a bimodal distribution. In addition, we restricted STDP to only effect synaptic 377 weights between excitatory neurons in our model. For inhibitory neurons, we used a normal distribution 378 of synaptic weights using a 10% variance. We scaled the synaptic weights of excitatory and inhibitory 379 connections from the starting distributions [0 1] to [0 4] for excitatory neurons and to [-14 0] for 380 inhibitory neurons. These scales were chosen to correspond to excitatory and inhibitory postsynaptic 381 potentials recorded in vivo (70). 382 For each network, we constructed the network topology and assigned synaptic weights between 383 connected neurons. We then allowed the network architecture to reweight with STDP until the firing 384 behavior of neurons to achieve activity that did not vary in firing rate or average oscillation rate by more 385 than 1% over a 5-minute simulation period. Across a range of connection architectures and synaptic 386 weights, we achieved a stable activity pattern for simulation times of at least 4 hours. For each condition 387 examined, we constructed and completed ten independent simulations, averaging the results from 388 these simulations into a single group. 389 We recorded four measures from each simulation: the neuron-activity rate, oscillation frequency, 390 oscillation peak magnitude, and oscillation width. Activity rate was calculated as the average firing rate 391 of neurons within the network over a five-minute period. To evaluate the occurrence of oscillations, we 392 recorded the number of action potentials occurring in a 10-millisecond sliding time window over the 393 simulation time ( Figure 6D). From these data, we selected times where peaks in activity occurred for the 394 network (peak prominence ≥ 1) and used these times to compute the oscillation frequency as the 395 number of oscillations per second. At each oscillation, we defined the oscillation magnitude as the 396 maximum number of spikes within the sliding window and the width as the full width at half peak 397 intensity. 398 oscillations per second, the peak number of spikes within our oscillation as a fraction of our network, 406 and the full-width-half-maxium (FWHM). (F) After network weights settled with STDP, we injured 407 networks at damage levels from 5% to 95%. After recording activity metrics immediately after injury, the 408 networks restructured connectivity weights, again according to STDP. We then reassessed network 409 function. 410

411
We further analyzed each network oscillation to identify preferred spike timings for neurons. For each 412 oscillation, we identified an interval that was 1.4 times the size of the oscillation width centered around 413 the peak. We then split the interval into uniform deciles and assessed the likelihood for each neuron to 414 fire within each time interval. 415 Damaging the neural network 416 Random neuron injury: To mimic traumatic injury, neurons were randomly removed from the network 417 after initial network reweighting with STDP. To maintain the initial excitatory-inhibitory balance of the 418 circuit, 4 excitatory neurons were removed for every 1 inhibitory neuron. To assess the immediate 419 effects of damage, we ran simulations without adjusting connectivity weights for 5 minutes, recording 420 both the average firing rate and oscillation parameters over this time period. Next, we allowed networks 421 to remodel with STDP with simulation times long enough to allow the firing rates and oscillation 422 behavior to settle (4 hours). We then reassessed our metrics with a stable connectivity for an additional 423 5 minutes ( Figure 6E). 424 Activity-based excitatory neuron injury: A second type of deletion scheme used the firing rate of 425 individual neurons as the selection criteria for deletion. Once a given network restructured with STDP, 426 we rank ordered the firing rate of each neuron within the network and deleted either the neurons with 427 the lowest firing rate (LFR) or, alternatively, the neurons that displayed the highest firing rate (HFR). 428 Because the focus of this removal strategy was to determine the effects of activity-based deletion, we 429 opted to only remove excitatory neurons. Similar to random neuron injury, we assessed the immediate 430 effects of neuron removal with a static connection weights for 5 minutes and then allowed the network 431 to resettle for 4 hours. 432

Statistical testing
To compare average activity and oscillation parameters between injury levels and baseline, we used 434 one-way analysis of variance (ANOVA) with Tukey-Kramer post hoc test. Comparisons between injured 435 and damaged networks at each injury level used t-tests with Bonferroni correction for multiple 436 comparison. 437