Computational model of neurodegenerative damage in the olfactory bulb: role of bulb and damage geometry, damage site, and identification of a potential non-invasive marker of early Parkinson’s disease onset

The olfactory bulb (OB) is one of the first regions of the brain affected by Parkinson’s disease (PD) as measured by both dysfunction and presence of α-synuclein aggregation. Better understanding of how PD affects OB function could lead to earlier diagnosis and potential treatment. By simulating damage to the OB in a computational model, it may be possible to identify regions of interest or markers of early disease. We modified a simple rate-based computational model of the olfactory bulb and simulated damage to various components of the network. This was done for several configurations of the network, at different sizes and with 1D and 2D connectivity structures. We found that, in almost every case, activity of 2D networks were more robust to damage than 1D networks, leading us to conclude that a connection scheme of at least 2D is vital to computational modeling of the OB. We also found that certain types of damage (namely, seeded damage to the granule cell layer and to the synapses between mitral and granule cells) resulted in a peak of the oscillatory power of the network as a function of damage. This result is testable experimentally and bears further investigation utilizing more sophisticated computational models. If proven accurate, this rise in oscillatory power in the OB has the potential to be an early marker of PD. Author summary One of the first symptoms of Parkinson’s disease is the degradation of the sense of smell. The olfactory bulb is the first region of the brain to process odor information and is affected by Parkinson’s disease at early stages. We simulated neural activity in a computational model of the olfactory bulb in the presence of damage and compared it to simulations of undamaged activity. We found that 2D model networks were more robust to damage than their 1D counterparts. We also found that 2D networks displayed increased oscillatory activity when damage was applied to certain parts of the network. This last result, if proven correct, would potentially be a marker of early-stage Parkinson’s disease, and if so, could aid in early diagnosis and treatment of the disease.

Introduction 1 excitatory signals to olfactory cortical regions. They also excite the interneurons of the 24 OB, called granule cells, which in turn inhibit the mitral cells through graded 25 dendrodentritic interactions [20]. This mix of inhibitory and excitatory interactions 26 drives oscillatory behavior in the OB network in the gamma band   [21]. 27 While the odor identity is first encoded in the OB, exactly how it is represented is 28 an ongoing question for which there are various theories, mainly revolving around 29 combinatorics of glomerular or mitral cell activity [22][23][24][25]. Oscillations in neural activity 30 in the bulb may also play a part in encoding odor identity, and are likely important to 31 odor recognition or information transfer, or both [21]. The precise manner in which PD 32 impacts OB function is still a matter of investigation [26,27], although misfolded 33 α-synuclein is present in the mitral, external tufted, granule, and periglomerular cells of 34 PD patients [5,28,29]. 35 It remains unclear how α-synuclein aggregates lead to neurodegenerative damage, 36 although a likely candidate for early cellular damage is disruption of transport along the 37 microtubule bundles of the axon with subsequent presynaptic degradation [30][31][32][33][34][35][36]. 38 There have been few investigations of the impact of damage that might arise from 39 neurodegeneration or lesions on neural networks, and most have examined the effect on 40 the classical Hopfield model of associative memory, in particular how memory capacity 41 is affected [37][38][39]. One paper examined how compensation could retain capacity in 42 such an associative memory model [40]. Importantly, we are not aware of any works on 43 biologically plausible sensory neural networks that examine the impact of 44 neurodegenerative damage with an intention towards guiding early diagnosis of the 45 diseases. 46 We investigate here the impact of neurodegenerative damage on various constituents 47 of a computational model of the OB network. Many excellent and insightful 48 computational models of the OB exist [41], focusing on various aspects of the olfactory 49 system, such as generating oscillatory behavior [42][43][44], glomerular layer 50 computations [45], and odor computations and representation [24,46]. As a simple 51 model of the olfactory bulb network-level activity that captures the oscillatory activity 52 in the gamma band [47], we chose to work with the Li-Hopfield model [42]. It simplifies 53 the OB to a network (with periodic boundary conditions) of mitral cells and granule 54 cells (the primary neuronal interaction in the bulb) (see Fig 1A). These cell types are 55 also of particular interest because some studies suggest that mitral and granule cells 56 August 1, 2020 3/30 may be especially vulnerable to damage in PD [9,29]. 57 The original Li-Hopfield work generated gamma band oscillations using a 58 one-dimensional network with ten mitral and ten granule cells, with an assertion that 59 similar behavior can be found with a two-dimensional network, and an implication that 60 the behavior was similar for larger networks. In fact, as we show here, the normal 61 behavior of a two-dimensional network is generally more robust against the model 62 damage as measured by how the odor representation diverges from the undamaged case, 63 and also shows, for some types of damage, a significant enhancement of gamma band 64 oscillatory power with increased damage. Because the real mitral-granule layers in the 65 olfactory bulb are two-dimensional, this represents an important difference between the 66 one-dimensional and two-dimensional versions of the model. In particular, we find that 67 damage at the synaptic interface or internal damage to the granule cells lead to 68 enhanced gamma band oscillatory power and a change in odor-evoked response pattern 69 at relatively larger levels of damage. In contrast, damage to the odor inputs or internal 70 damage to the mitral cell inputs lead to reduced oscillatory power with damage but 71 shows a similar modification of the odor-evoked response with damage in 72 two-dimensions. In one-dimension, no matter what the region of damage, typically the 73 gamma band oscillatory power drops with damage, and the odor representation diverges 74 from the undamaged one more quickly than in two-dimensions. respectively. The exact functions and values for the parameters can be found in S1 94 Table. 95 The connection matrices H 0 and W 0 dictate the structure of the network. The By taking a computational model of the olfactory bulb and perturbing the network, we 116 can make predictions regarding the impact of various types of damage. The simulations 117 are run in Python using Scipy's solve ivp function [48], and all code for the simulations 118 can be found on the Github repository (see Supporting Information). Damage to the network is measured by δ, the fraction of weight removed. For the weight matrices, W 0 120 for example, where W 0 is the undamaged matrix and W Damaged is the damaged matrix. In a given 122 trial, damage is delivered to one of the following components of the model: implemented by multiplying the right-hand-side of the differential equation (except for 130 the leak term) for the given cell by some fraction less than one, (1 − δ i ). For example, 131 would be damage delivered to the i th mitral cell unit. In this case, δ is calculated as where N is the number of mitral cells. 133 We ramp up the damage to the selected part of the network, run the network at that 134 damage level, and compare the activity to the activity of the undamaged network.  For FD, the damage is delivered to every element of the selected component equally. 138 This amounts to simply scaling the chosen quantity uniformly. For example, if FD was 139 applied to H 0 , each element of H 0 would be reduced by the same fraction of its original 140 value on each damage step. This continues until the matrix is reduced to zero (see Fig 3 141 A).

142
For CD, the damage is delivered to a specific element (or column if the network 143 quantity is a matrix), ramped up until that element is reduced to zero, and then that 144 procedure is repeated on an adjacent element until the maximum damage level is 145 reached. For example, if CD was delivered to H 0 , damage would be delivered 146 incrementally to a single column (i.e., the postsynaptic strength of a single granule cell) 147 until it was reduced to zero. The same process would then begin on the column to the 148 right (see Fig 3 B). We continue this process until 50% of the columns are removed.

149
SD is a combination of FD and CD. Damage is first delivered to a single element (or 150 column), and on the next damage step, damage is delivered to that element again, as 151 well as to neighboring elements. For example, if SD was enacted on H 0 in a 1D network, 152 it would begin on one column, say column 6. On the subsequent damage step, the 153 damage would be delivered to columns 5, 6, and 7. This spreading continues with each 154 damage level until the matrix is reduced to zero (Fig 3 C). Characterizing Network Activity 156 We characterize the damaged network's activity by comparing its activity pattern, 157 average activity level, and average oscillatory power to that of the undamaged network. 158 The activity pattern is measured by low-pass filtering the cell activity and averaging over time for each mitral cell [42]. This yields a response vector R, where each entry reflects the average activity level for a single mitral cell unit. We then compare the response vector for the damaged network, R D , with the response vector for the undamaged network, R 0 , by calculating the distance between activity patterns, D P [42]: This normalizes the maximum distance between patterns to D P = 1, with identical 159 patterns giving D P = 0. If D P increases with damage level, that indicates that the 160 damaged network's activity pattern for a given odor input is diverging from that of the 161 undamaged network.

162
The average activity level, R, is calculated as a root-mean-square of the response vector R. We measure the difference between average activity levels, D A , by [42]: which has a max of |D A | = 1 and a min of |D A | = 0. We retain the sign in practice 163 because it indicates whether the average activity level of the damaged network is greater 164 or less than the null case. 165 We calculate the average oscillatory power, P avg , by first high-pass filtering the 166 mitral cell activity above 15 Hz to ignore theta band (2-12 Hz) activity, which was 167 beyond the scope of the present study. We next calculate the power spectrum (P(f )) 168 for each mitral cell from 125 ms to 250 ms using Scipy's periodogram function [48]. This 169 time window captures the oscillatory behavior during the most active part of the cycle 170 (see 1C) while ignoring the spurious signals that can arise at higher levels of damage 171 that are not actually due to gamma band oscillatory activity (see S10 Fig). We then 172 integrate the power spectrum over all frequencies, f , for each cell (see S10 Fig) and 173 average over the mitral cell population (N ) to get P avg ., Each quantity is averaged over five trials. For CD and SD, we then repeat the trial 175 using each cell as the starting point and average the values again over all starting cells. 176

177
Damage to W 0 , H 0 , and GCL 178 CD caused network activity to diverge from the undamaged case at the lowest δ, 179 followed in impact as a function of δ by SD and FD respectively, as shown in  This was consistent across network sizes and structures for damage delivered to W 0 , H 0 , 181 and GCL (see S1 Fig -S6 Fig). The 2D networks were more robust than their 1D counterparts in almost every case, 183 as measured by how quickly their activity diverged from the undamaged case (Fig 5). 184 This is especially apparent in SD delivered to W 0 , H 0 , and GCL, but can also be seen in 185 CD and FD for those neuronal regions (S1 Fig -S6 Fig). The effect of damage on average oscillatory power also depended on the damage 187 scheme and network structure. FD and SD in 2D networks resulted in increases in P avg 188 at intermediate δ (Fig 6A), but CD rarely showed an increase in P avg (S8 Fig). While 189 this rise in P avg was ubiquitous among 2D networks with FD and SD, P avg decreased 190 monotonically for most 1D cases (see Fig 6B).  Distance in activity pattern for SD delivered to W 0 , H 0 , and GCL in the 1D and 2D 100 cell networks plotted against damage as a fraction of total synaptic weight removed (damage level, δ). In each case, D P increases to its maximum value most rapidly for the 1D network, indicating a greater divergence from the activity pattern in the undamaged case at lower levels of damage. B: Distance in average activity level for SD delivered to W 0 , H 0 , and GCL in the 1D and 2D 100 cell networks. In each case again, the 2D network activity is more robust to damage.

192
Trials with damage delivered to MCL and to OI showed qualitatively similar results.

193
Like the W 0 , H 0 , and GCL trials, network activity patterns were most sensitive to CD, 194 followed respectively by SD and FD (Fig 7A). However, the average activity level 195 changed at a similar rate for all types of damage propagation, though usually most 196 slowly for CD (Fig 7B). Additionally, P avg decreased in every case at similar rates 197 (Fig 7C).

199
The work described here illustrates the importance of the 2D network structure of the 200 OB in terms of robustness against neurodegenerative damage (relative to 1D models), 201 Fig 6. Effect on Oscillatory Power. A: Average oscillatory power (P avg ) for FD, CD, and SD delivered to H 0 in the 2D 100 cell network plotted against damage as a fraction of total synaptic weight removed (damage level, δ). B: P avg for SD delivered to H 0 in the 1D and 2D 100 cell networks. FD and SD to H 0 result in a rise in P avg for the 2D network, while CD to the 2D network and any kind of damage to the 1D network do not. Effects of Damage to Mitral Cells or to Olfactory Input A: Distance in activity pattern for FD, CD and SD delivered to MCL and OI in the 2D 100 cell network plotted against damage as a fraction of total synaptic weight removed (damage level, δ). B: Distance in average activity level for FD, CD and SD delivered to MCL and OI in the 2D 100 cell network. As opposed to damage to W 0 , H 0 , or GCL, average activity level very nearly decreases monotonically for all cases, with the exception of CD. C: Average ocillatory power for FD, CD, and SD delivered to MCL in the 2D 100 cell network.
as well as in modification of the power spectrum with damage. Computational models 202 of the olfactory bulb have synaptic connection schemes that are mainly either 2D or 3D 203 (for example, [44,46,49]), or random (for example, [43,50]) in structure. We confirm 204 that for applicability to disease models, accounting for the geometric structure of the 205 bulb is important, as even in the simplistic model implemented here there was a sizeable 206 difference in response to damage based on the dimensionality and geometry of the 207 connection scheme.

208
The structure of the synaptic connections also has implications for the spread of 209 pathology. We propagated damage here in three ways, FD, CD, or SD. FD and CD 210 represent the two extremes of damage propagation (FD completely delocalized to every 211 cell, and CD completely localized to a single cell before spreading). Misfolded 212 α-synuclein spreads from cell to cell before cell death [51,52] in what many believe is a 213 prion-like manner [7,13,53], although the precise mechanism and the level of damage in 214 the donor cell before transmission is still under investigation [54,55]. This means the 215 damage would spread differently for a 1D connection scheme compared to 2D. It also 216 makes the SD method the most relevant form of propagation presented here, though FD 217 and CD are still important for comparison.

218
The results presented above imply that localized damage may be more severe than 219 dispersed damage at similar levels as measured as a fraction of synaptic strength, 220 though some experimental results do not support this [56]. This may be a failing of the 221 model, and taking the plasticity of the bulb into account would likely alter these results. 222 The work here also illustrates a differentiation between two general groups of 223 damage types. Assuming SD is the most realistic propagation strategy, damage to H 0 , 224 W 0 , and GCL all result in increased network activity and a peak with damage in 225 gamma band oscillatory power, while damage to MCL and OI do not. This gives 226 specific differentiating factors that can be measured.

227
One of the most salient results presented here is the peak in P avg for SD to W 0 , H 0 , 228 and GCL for 2D networks. This provides a testable marker that can help to identify 229 potential network components damaged in PD. In fact, Kulkarni et al. [27] recently 230 recorded local field potentials in mice olfactory bulbs after injection of α-synuclein 231 pre-formed fibrils. They saw a rise in oscillatory power during odor presentation after in 232 incubation period ranging from 1-3 months post-injection. This increase in oscillatory 233 power was in the beta band only (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30), while this network models activity in the 234 gamma band , so comparison is limited. However, beta oscillations are likely 235 mediated by the same synapses in the OB as gamma [57,58], and we believe these 236 results to be related. Taken together, they imply PD may damage the synapses between 237 the MCL and the GCL (in this model, W 0 and H 0 ), although generalized damage to the 238 GCL shows similar results here. Studies have shown that α-synuclein plays a role in 239 synaptic transmission [59,60], and that α-synuclein oligomers and aggregates disrupt 240 synaptic function [59,[61][62][63]. Additionally, there is evidence that neurodegeneration in 241 PD may be due to synapse dysfunction rather than just cell death [64]. Therefore, we 242 believe the W 0 and H 0 results to be the most interesting. They are particularly 243 intriguing since a recent study has shown that it is possible to measure such an increase 244 in oscillatory power without an invasive probe to the brain [65], so this can be regarded 245 potentially as a unique early onset marker of the disease.

246
Opportunities exist for experiments to investigate the results and claims presented 247 here. More acute measurements of olfactory bulb oscillatory activity early on in 248 pre-formed fibril seeding before severe cell damage has a chance to occur could help 249 illuminate whether effects are due to synaptic dysfunction or cell death and could lead 250 to more directed studies in the future.

251
The model employed here is simplistic. Repeated trials using more biophysically 252 relevant models of the olfactory bulb should be performed. For example, future work 253 should implement a model with activity in both the gamma and the beta bands, such as 254 in models by Osinski & Kay [50] and David [66], requiring more sophisticated 255 centrifugal inputs [57]. Comparing responses to damage across models can uncover 256 commonalities and make reasonable predictions about actual pathology and behavior in 257

Supporting information
All code for the simulations can be found at https://github.com/jkberry07/OB PD Model. S1 Table contains the equations and parameters used in the equation. For more information about the parameters, and for information about the noise added to the system, see the original paper by Li and Hopfield [42], as well as Li's dissertation [67]. Supplemental tables and figures are located after the references. S1 Table. List of parameters and functions. Cell activity (g x ) at various levels of seeded damage to W0 in the 2D 100 unit network, with the associated power spectrum plotted next to each. The inhale peak is at 205 ms. The oscillation amplitude can be seen to grow from δ = 0 to δ = 0.8448, which is reflected clearly in the power spectra. Note that the vertical scale is larger in the last cell activity plot. The sharp fall in activity from about 225 ms to 300 ms can cause an increase in P avg , but clearly is not actually an example of gamma band oscillatory behavior and is the reason for the shortened time window used for August 1, 2020 28/30 calculating the power spectra for P avg . As stated in the text, the power spectrum was calculated from the cell activity high-pass filtered above 15 Hz from 125 ms to 250 ms. Some power density was present above 100 Hz, but power density below 100 Hz dominated the contribution to P avg .

S11 Fig. Evolution of Cell Activity in Flat Damage to H0 in the 1D 100 Unit Network
Cell activity (g x ) at various levels of flat damage to H 0 in the 1D 100 unit network, with the associated power spectrum plotted next to each. The inhale peak is at 205 ms.