Morphological variability may limit single-cell specificity to electric field stimulation

Non-invasive brain stimulation techniques are widely used for manipulating the behaviour of neuronal circuits and the excitability of the neurons therein. While the usage of these techniques is widely studied at the meso- and macroscopic scales, less is known about the specificity of such approaches at the level of individual cells. Here we use models based on the morphologies of real pyramidal and parvalbumin neurons from mouse primary visual cortex created by the Allen Institute for Brain Science to explore the variability and evoked response susceptibility of different morphologies to uniform electric fields. We devised a range of metrics quantifying various aspects of cellular morphology, ranging from whole cell attributes to net compartment length, branching, diameter to orientation. In supporting layer- and cell-type specific responses, none of these physical traits passed statistical significance tests. While electric fields can modulate somatic, dendritic and axonal compartments reliably and subtype-specific responses could be observed, the specificity of such stimuli was blurred by the variability in cellular morphology. These null results suggest that morphology alone may not account for the reported subtype specificity of brain stimulation paradigms, and question the extent to which such techniques may be used to probe and control neural circuitry. Author summary Over the last several decades there has been increased interest in the efficacy of non-invasive brain stimulation, particularly in determining the limits of specificity of such techniques. Despite this growing area of research, much remains unknown about the interactions of non-invasive techniques with neurons at the single-cell level, notably the importance of morphology to these interactions. We make use of detailed single-neuron models and simulate them in a uniform electric field and demonstrate that the high variability in neuron morphologies may limit how specifically single neurons can be targeted non-invasively. We confirmed this for neuron morphology characteristics at macro- and meso- scales and at varied orientations. Our work suggests that previously reported subtype specificities in non-invasive frameworks are not accounted for by considering only morphological factors.

Non-invasive brain stimulation (NIBS) paradigms are used for a variety of purposes, 2 including attempted treatment of several psychiatric conditions (e.g. major depressive 3 disorder [1,2], schizophrenia [3], bipolar disorder [4]), and investigating neural 4 oscillation properties and controllability [5][6][7]. Despite reported clinical efficacy, many 5 questions remain about the mechanisms involved [8] as well as how to optimize their use. 6 The high variability in brain response to NIBS protocols represents a major limitation 7 in these efforts. Such variability has been attributed to genetic and neurophysiological 8 factors, but evidence suggests brain state during stimulation may further impact the 9 response [9][10][11]. Many contributing factors to response variability are immutable (e.g. 10 subject age, skull thickness, etc), but others, such as the placement of the apparatus (e.g. 11 electrodes or coil), are malleable [10,12]. The immense combination of free parameters 12 in these protocols greatly complicates parsing which parameters are crucial. This has 13 led to the development of a variety of computational models in efforts to improve the 14 efficacy of experimental approaches and our overall understanding of NIBS [7,13]. 15 Despite widespread utility, comparatively little is known about the effects of NIBS at 16 the level of single neurons [8]. It is known experimentally, although not widely 17 investigated, that NIBS techniques can modulate the activity of individual 18 cells [12,14,15]. Along with this awareness has come an acute interest in parsing the 19 limits of specificity of these techniques. That is, the extent to which neuron traits (e.g. 20 subtype, morphology) govern the ability of a single NIBS protocol to acutely activate or 21 suppress the response of neurons sharing similar characteristics. Preliminary efforts 22 have been made to elucidate a relationship between response variability and neuron 23 morphology via intermediate morphological models (that is, models with a limited set of 24 stereotyped, structural components) with nominal success. One such simplified model 25 examined the relationship between many physical traits (e.g., branching, compartment 26 widths and lengths, etc.) of a basic neuron and found, on a trait-by-trait basis, they 27 alter the strength of electric field (E-field) necessary for the neuron to fire [16]. 28 Moreover, even for the most simplified case of a straight-cable neuron model [17], 29 changing a single property, such as the length of the cable, was found to influence the 30 response. Collectively this supports the consensus that physical morphology acts as an 31 avenue for specificity within NIBS frameworks. 32 Experimental work in alert non-human primates, at field strengths typically applied 33 to humans, has shown transcranial electrical stimulation to reliably control the timing, 34 but not firing rate, of spikes in individual neurons [14]. However, that experimental 35 work found no notable difference in phase locking between pyramidal (PC) and 36 parvalbumin (PV) cells, a result at odds with what has historically been expected of 37 individual neurons, which are long theorized to have subtype-and layer-specific 38 differences in NIBS response [7,18,19]. In light of this, it may then be purported that 39 other neuron traits (e.g. morphology) that vary widely across neuron populations may 40 act as the limit to specificity. 41 Here, we put these results to the test using detailed morphological 42 multi-compartment models from real mouse neurons to computationally study the 43 effects of NIBS on evoked responses of neurons of different types (PV, PC) and 44 populating various cortical layers via the application of a uniform E-field. We 45 specifically examined which, if any, physical characteristics might be involved in 46 generating cell-type and/or layer-specific responses. To do so, we developed metrics 47 quantifying various aspects of cellular morphology, both at the level of whole cells and 48 the compartments they are modelled from, ranging from cell volume to compartment 49 lengths, diameters, number of branches, cellular orientation as well as interneuron 50 myelination [20][21][22] The models used here were created by the Allen Institute for Brain Science based on the 66 morphologies and biophysical properties of real neurons found in the primary visual 67 cortex of mice [20,21]. That is, these computational mathematical models of single    [20,25,26]. Subsequent data analysis was performed using the NumPy, 92 Matplotlib, Seaborn, SciPy, and Pandas Python libraries [27][28][29][30][31]. 93 In the NEURON environment, the membrane potentials of the models are calculated 94 using the cable equation: where I net and V are the net current (ionic and injected) and membrane potential, 96 respectively. This is then approximated to its spatially discretized form such that the 97 neuron is reduced to a set of connected compartments, and Eq (1) becomes a family of 98 equations, where c j is the membrane capacitance of compartment j, i inj,j are the injected 100 currents, and i ion,j includes all currents through ionic channel conductances. The contains all of these compartments as segments of the section [25,32]. This is done for 108 computational efficiency, as Eq (2) may then be re-formulated as, where ∆x and d are the compartment length and diameter, respectively. If these are 110 the same between compartments, as is the case for a section, then R a ∆x/π(d/2) 2 is the 111 axial resistance, and C m πd∆x is the compartment capacitance.

112
The electric field is applied to all compartments in a cell using NEURON's 113 extracellular function [25]. This function adds the field in millivolts such that it is 114 connected in series with the conductance of the last extracellular layer of the model. For 115 these simulations, we use a uniform electric field in the positive z-direction, such that 116 each compartment receives field E z = E cosθ + ξ, where θ is the angle between the 117 middle of a given compartment and the z-axis, and ξ is noise. Ten trials each of 118 subthreshold field strengths E = {−50, −30, −10, 0, 10, 30, 50} mV were applied, in line 119 with ranges that have been used for similar models [16,18]. The Gaussian white noise, ξ, 120 was added with the field applied to each compartment and is different between 121 compartments and across trials. The cell response was read out from the somatic 122 compartment as a membrane potential and the average response across trials was taken 123 to quantify the specificity of the cell.

124
Characterizing neural morphology 125 In the present framework, the Cartesian coordinates for all compartments in the models 126 are known. This allows for both whole-cell and compartmental characterizations of the 127 morphologies to be made. For macroscopic measures, that is characterizations that the number of compartments in model j. For the z -direction length (referred to 133 hereafter as length), the value L is calculated as simply the difference between the 134 maximum and minimum z-direction coordinates. Finally, the elliptical volume is 135 calculated by taking the length, and the equivalent maximum distances in the x -and y-136 directions, halving each of them and utilizing the formula for the volume of an ellipsoid, 137 V = π/6(x max y max z max ).

138
Additionally, at the level of whole cells, we sought to investigate the effects of 139 orientation. To do this, the models were rotated about the y-axis using a rotational work on these neuron types [22].

168
Susceptibility 169 The curves of membrane potential resulting from the application of various applied field 170 strengths as recorded at the somatic compartment in the upright (or 0 • ) orientation 171 serve as a representation of the susceptibility of a neuron to the applied field. Here, we 172 measure this property using the slope of the membrane potential curves: where ∆E is the difference between the E-field strengths for the cases (here equal to 174 100) and ⟨V ⟩ i are the average membrane potentials at field strengths i = {−50, +50}.

175
As ∆E = −100 here, the resulting values of S are also negative. calculates the Pearson correlation coefficient, R, and its corresponding p-value using a 181 Wald Test with t-distribution of the test statistic [30]. These correlations then allow us 182 to discern which morphological trait (if any) may contribute to the specificity of the 183 NIBS protocol.

185
Characterizing neuron morphology variability been demonstrated that, in isolation, each of these features may affect the field strength 215 required to evoke depolarization [16]. For these models, the compartment lengths show 216 more variability across the layers and subtypes in the axonal compartments than in the 217 dendritic ones (Fig 2A). Further, the dendritic compartments are on average longer 218 than the axonal compartments, however, there are exceptions such as in L5 PC models 219 where axonal and dendritic compartments were around the same length.

220
By contrast, the compartment diameters are relatively similar between the dendritic 221 and axonal compartments (Fig 2B). The diameters observed between subtypes and 222 layers are also, on average, relatively similar to each other. Branching is, generally, 223 observed to be much more prevalent in the dendritic compartments than the axonal 224 ones (Fig 2C). However, there are exceptions to this, such as in the L5 and L23 PV  Subtype-and layer-specific evoked responses to uniform electric 229 fields 230 We first examined the evoked response of these models across layers and subtypes, 231 independently of their morphological variability. To do this, we apply a uniform electric 232 field to the neuron models (Fig 3A). The sign and strength of the field will influence hence the relative evoked response may be quantified by averaging membrane potential 235 fluctuations across independent trials (see Fig 3B and Methods section).

236
While not statistically significant, differences in evoked responses could be observed 237 between the subtypes and layers. Pooling the neurons by layer (neglecting the subtypes) 238 there is no significant difference between the somatic responses observed. The lowest 239 coefficient of variation of these membrane potentials, ⟨CV V ⟩, is observed for the L23 240 models, with the highest being the L4 models (Fig 3C). By comparison, ignoring the 241 layers and pooling the neuron models based on their subtypes, there are visible 242 differences between the depolarization responses of PC and PV cells (Fig 3D). Of note, 243 the differences observed between subtypes manifest differently in the cases where the 244 PV neuron models are myelinated or not. For unmyelinated PV neurons, the same field 245 strength leads to less depolarization than is observed in the PC models at the same field 246 strength. However, in the myelinated PV models, the response to the same field 247 strength eventually catches up with the PC models in terms of depolarization.

248
Interestingly, the myelinated PV models also hyperpolarize more aggressively than 249 either the PC or unmyelinated PV models (Fig 3D). The lack of statistically significant 250 differences observed in these poolings suggests limitations to the extent of specificity 251 attainable via NIBS techniques for circuit control. The average CV V for the pooled 252 subtypes suggests that the PC models have more variation in their responses. Based on 253 the severity of this curve compared to the PV models and the similarly shaped curves in 254 the same plot pooled by layers, it is likely that the majority of this variation comes from 255 the L4 PC models.   (Fig 4). That is, while rotating the neuron by 180 o , or even 90 o , will 265 result (for this neuron) in a stronger depolarization however, the net response to the 266 applied field is largely conserved at the soma, and do not differ significantly from the 267 somatic readings of the other models (shown in grey). However, it does have a marked 268 effect on the response of a given dendritic or axonal compartment (see Fig 4B middle   269 and right). This suggests that, for neurons in these layers, the attainable specificity 270 through NIBS techniques is insufficient to control neurons with a specific orientation. Previous NIBS studies have demonstrated that single neurons are affected and entrained 294 by low field strengths [14,17,[33][34][35]. The extent of specificity NIBS holds over single 295 neurons is less clear and, despite some successes that found no difference in the 296 entrainment of excitatory and inhibitory neurons [14,15], is largely experimentally  2), the slope of the average somatic membrane potential curves is shown as a function of D) average compartment length, E) average compartment diameters, and F) number of branches. All slopes are calculated from response in the upright (0 • ) orientation. Schematics to the right of each plot correspond to the quantity on their x-axis. The Pearson correlation coefficients, R, their significance, p, and the line of best fit are found using SciPy's linear regression package [30]. variability and physical morphology. These results align with experimental work 300 revealing that highly complex, variable morphologies had no effect on the conservation 301 of physiological waveforms and circuit functions of neurons [36]. Similarly, recent 302 computational work in L5 pyramidal cells has shown morphological variance to be 303 insufficient to reproduce electrical variability [37].

304
Among the morphological diversities we considered was myelination, which influences 305 neuron response to NIBS [38][39][40]. The minor differences in the PV, but not PC, models 306 in the myelinated and unmyelinated versions (see Fig 3D) suggest that response 307 variability due to myelination is non-uniform and, consistent with a recent work [33], 308 requires >15% coverage to manifest. This may contribute to earlier computational 309 works, without myelination, that found layer-specific differences in the E-field strength 310 required to depolarize neurons in L23 and L56 [18] whereas recent studies involving 311 myelination have shown minimal difference required to evoke firing from L23, L4, and 312 L5 neurons [33,41]. As the effect of myelination is related to its amount, in brain 313 regions with less myelin, or shorter axons, inhibitory neurons may respond more 314 similarly to their unmyelinated equivalents (see Fig 3D), and hence be less depolarized 315 than excitatory neurons. Indeed, such effects may explain differences in our results and 316 other myelinated models that predicted specific activation of cells across subtypes and 317 layers. Where those studies use morphologies from multiple regions of the brain [19], 318 they may indicate that morphologies across brain regions differ more substantially than 319 those found in the same region, and are hence more specifically selected. Such a 320 distinction may also explain discrepancies between those works and the lack of 321 significant difference found between subtype responses in experimental protocols [14].

322
In addition to layer-and subtype-based neuron groupings, we sought to explore the 323 influence of physical morphology on neuron response (see Fig 5). Compartment-level 324 studies of simplified neuron models previously found individually, varied physical traits 325 which influence the field strength required to evoke firing [16]. Based on those findings, 326 one could infer that an ideal neuron (that is, one with the lowest required field strength 327 required to fire) would have long, thick dendrites with many branches that disperse 328 minimally from the primary axis the applied field lays on. Moreover, an ideal neuron 329 would have long, narrow axons with minimal branching that disperses widely from the 330 field axis. For the models considered here, there is no class of neurons that fits this ideal 331 (e.g. long but thicker axons; see Fig 2). In line with this, our simulations did not 332 demonstrate any singular physical trait to be the driving source response specificity to 333 the same stimuli between neurons. This suggests the difference between individual 334 neuron models and their response is likely the result of a convolution of the physical 335 properties rather than individual traits determining the outcome.

336
For computational efficacy, many models of NIBS protocols make use of mean-field 337 type approaches, whether for the whole model population [42] or for subtypes within 338 the population [7]. If groupings of neurons could be derived on the basis of physical 339 morphologies, this would imply the need for adaptations to these models to include 340 more extensive morphological features. In contrast to this, here we have found no 341 significant relationship between susceptibility and morphology ( Fig 5); this suggests 342 support for the assumptions made in designing such computational models to not 343 include specific morphology in their frameworks.

345
Our morphologies come from neurons found in the primary visual cortex of mice; 346 importantly, this limits the extrapolation of the results of these simulations to any 347 expectations for human cells, as it has been shown that a number of neuron properties 348 do not scale from non-human to human neurons [43,44]. In line with this, recent results 349 showed human L5 PC neurons have unexpected biophysical properties for their size, 350 despite the composition and allometry of the human cortex scaling relative to the nine 351 other species tested [44]. Further, as these models are simulated in isolation, they 352 disregard any potential influence from conductive tissue, distance-to-skull or other cell 353 interactions that may bolster response through biophysical mechanisms, such as 354 ephaptic coupling [13,34,35,45]. White noise is added to account for some of these 355 effects, however, the net efficacy may still be impacted. Additional model-specific 356 limitations can be found in the Allen Institute technical papers and 357 documentation [20,21].

358
The neuron models were predominantly simulated in one orientation in free space.

359
However, the overall effects of orientation were considered by rotating the neuron 360 models 0 • , 90 • and 180 • about the y-axis and recording from the somatic, axonal and 361 dendritic compartments in all three orientations (see Methods and Fig 4). The observed 362 effect of this change is relatively small at the somatic compartment, although the overall 363 polarization of the model changes. NIBS studies of biophysical models found this 364 minimal effect results from the axon branching and myelination, which drastically 365 reduces the effects of orientation on activation threshold [33]. The layers considered may 366 also reduce the effect of orientation, as L1 and L6 have previously shown more 367 variability in preferential orientation [8,41]. The average membrane potential in 368 response to orientation changes is more attenuated at the soma than the dendrites and 369 axons, which are two to three times more susceptible to polarization at their terminals 370 than somas [46]. 371 Importantly, for the model neurons, there is no distance component to the applied 372 field as it is applied as directly to the surface of the neuron (see Methods). For real 373 neurons receiving input from outside the skull and encircling brain tissue, the minimum 374 field strength required to facilitate depolarization may be greater than in these models 375 where the field is applied uniformly. Alternatively, having connections to other neurons 376 to bolster the response may also allow in vivo neurons that would not individually react 377 at a lower field strength to be pushed toward depolarization.

379
Non-invasive brain stimulation techniques are widely used for manipulating neuronal 380 circuit behaviour and excitability. While the usage of these techniques is widely studied 381 at the compartmental and whole-cell scales, less is known about the specificity of such 382 approaches at the level of individual cells. Using models based on real pyramidal and 383 parvalbumin neurons, we characterize their morphologies to demonstrate the observed 384 variability in physical traits limits the specificity attainable with electric field 385 stimulation. Our null results suggest that morphology alone may not account for the 386 reported subtype specificity of brain stimulation paradigms and that such techniques 387 may be limited in their specificity by morphological variability.