Inter-species h-current differences 1 influence resonant properties in a 2 novel human cortical layer 5 neuron 3 model 4

15 Most existing multi-compartment, mammalian neuron models are built from rodent data. However, 16 our increasing knowledge of differences between human and rodent neurons suggests that, to 17 understand the cellular basis of human brain function, we should build models from human data. 18 Here, we present the first full spiking, multi-compartment model of a human layer 5 cortical 19 pyramidal neuron. Model development balanced prioritizing current clamp data from the neuron 20 providing the model’s morphology, while also ensuring the model’s generalizability via preservation 21 of spiking properties observed in a secondary population of neurons, by “cycling” between these 22 data sets. The model was successfully validated against electrophysiological data not used in 23 model development, including experimentally observed subthreshold resonance characteristics. 24 Our model highlights kinetic differences in the h-current across species, with the unique 25 relationship between the model and experimental data allowing for a detailed investigation of the 26 relationship between the h-current and subthreshold resonance. 27

To explore this seeming inconsistency, a combination of computational and experimental 54 techniques are employed to create a novel human neuron model with a particular focus on the 55 h-current. The development of computational models of human neurons with high levels of 56 biophysical detail are more challenging than their rodent counterparts due to limited access to 57 tissue for experimental recordings. This challenge is exacerbated by the fact that to model how 58 a specific channel contributes to cellular dynamics, it is typically necessary to obtain a complete 59 data set (including whole-cell recordings in current and voltage clamp modes, pharmacological 60 manipulations, and 3D morphology) all in the same neuron. The increased access to rodent tissue 61 makes accounting for these concerns more feasible in the rodent setting, and explains why a 62 majority of the existing biophysically detailed neuron models are constrained by rodent data (Dong,  neuroscientists: in what settings is it appropriate to utilize rodent neuron models to glean insights 67 into the human brain, and when such approximations are undermined by inter-species differences, 68 can the functional role of these differences be identified? 69 We address these questions via a modeling framework that makes use of a detailed data 70 set obtained from a single human L5 neuron. We are motivated by the clear preponderance 71 of the h-current in human L5 neurons (Chameh et Ulrich, 2002), and towards developing human inspired neuronal models for brain simulators 75 (Einevoll et al., 2019). 76 Since it is clear that the characteristics of a given cell type are not fixed (Marder and Goaillard,77 2006), and moreover that this inherent variability amongst similarly classified cells could be func- 78 tionally important (Wilson, 2010), we develop a modeling approach that directly accounts for the 79 challenge posed to modelers by such "cell-to-cell variability". Our "cycling" model development 80 strategy primarily constrains the model using current clamp data and morphology from the same 81 neuron. In a second step, we ensure the model retains spiking characteristics exhibited by a popu-82 lation of secondary human cortical L5 pyramidal neurons; the process cycles between these two 83 steps to obtain an optimal model. The resulting multi-compartment, fully spiking human L5 neuron 84 model recapitulates the electrophysiological data from hyperpolarizing current clamp experiments 85 in the primary cell remarkably well, while also demonstrating repetitive and post-inhibitory rebound 86 spiking properties characteristic of human L5 pyramidal cells from the secondary data set (Chameh 87

Figure 1. Morphology and current clamp data obtained from the primary neuron. (A)
The morphology of the primary neuron was reconstructed using IMARIS software and imported into NEURON (which generated the plot shown here). (B) Current clamp recordings from the primary neuron in the presence of TTX that are the primary constraining data for model development.
carinic potassium channel (abbreviated Im). Note that the abbreviations used here are motivated 139 by the labeling used in the NEURON code for consistency. This provided the initial basis for our 140 model construction, with further details included in the Materials and Methods. 141 The "cycling" technique schematized in Figure 2 built upon this basis. In the first step, an 142 optimization algorithm was run to best "fit" the model's output with blocked sodium channels 143 to experimental data from current clamp recordings in the presence of TTX (see Figure 1). This 144 determined a majority of the conductances used in the model, as well as the passive properties 145 and the kinetics of the h-current. As the h-current is the primarily active inward current at hyperpo- kinetics of this channel type. 149 In the second step, after a best fit was achieved, we hand tuned the conductances involved in 150 action potential firing (sodium conductances and the K_Pst and SKv3_1 potassium conductances, 151 which were not altered in the preceding step), along with minor alterations to the dynamics of these 152 channels (see details in the Materials and Methods). The goal of this step was to ensure the spiking 153 behavior of our model cell was reasonable in comparison to the range of spiking properties, both of 154 repetitive and post-inhibitory rebound (PIR) firing, exhibited by secondary human L5 pyramidal cells 155 (summarized in Table 1 (Chameh et al., 2019)). We aimed to obtain these firing characteristics with 156 minimal potassium conductances, in order to minimize the error seen in Figure 3E: an extensive 157 exploration of the parameter space revealed that a "best fit" of this trace would enforce values of 158 the potassium conductances that would not permit action potential firing, motivating the hand 159 tuning of these values in search of a set of sodium and potassium conductance values that would 160 permit spiking while also minimizing this error. As the properties of these potassium channels 161 Figure 2. Diagram of the model development strategy. Hyperpolarizing current clamp data taken from the primary human L5 pyramidal cell was the primary constraint in determining model parameters. To ensure that the model exhibited repetitive and post-inhibitory rebound firing dynamics characteristic of human L5 pyramidal cells, data from secondary neurons, as well as best fit data from depolarizing current clamp experiments in the primary cell were used, and a "cycling" technique was developed in which conductances primarily active during spiking dynamics were fit separately by hand. The adjustments to the potassium conductances affect the current clamp fits, so these were re-run with the new values, hence the "cycle". affected the current clamp fits, it was then necessary to run the optimization algorithm of the first 162 step again with these new values, hence the "cycling". This cycling pattern continued until no further 163 improvement in the model, as determined via the quantitative error score from the optimization 164 process as well as the more qualitative matching of spiking properties, could be obtained (see the 165 Materials and Methods for further details). 166  in Figure 3A-E. The repetitive spiking behavior of the model in response to various driving currents is 169 shown in Figure 4A-C, and the capacity for PIR spiking is shown in Figure 4D; both of these protocols 170 are performed with active sodium channels. The repetitive firing frequency or latency to the first 171 PIR spike (depending upon whether the protocol is a depolarizing or hyperpolarizing current clamp, 172 respectively), is given in Table 1. Critically, the model closely matches all of the hyperpolarizing 173 current clamp data, indicating that the dynamics of the h-current within this voltage range were 174 accurately encapsulated by our model. While the error in the depolarizing current clamp recording 175 ( Figure 3E) is more noticeable, this was minimized via the process described above, and was the 176 best case while also ensuring reasonable repetitive spiking and PIR spiking behaviors (Figure 4 and   177   Table 1).
178 Figure 3. Model well fits data from hyperpolarizing current steps, in which the h-current is the primary active channel, while minimizing the error seen in a depolarizing current step. (A-D) Fits of current clamp data with -400 pA (A), -350 pA (B), -300 pA (C) and -50 pA (D) current steps with TTX. (E) Fit of current clamp data with a depolarizing current step of 100 pA with TTX. All four hyperpolarized current steps are fit with great accuracy, with a focus on the initial "sag" and post-inhibitory "rebound" that are driven by the activity of the h-current. While the charging and discharging portion of the depolarizing current trace is well fit, the amplitude of the response is less accurate; however, this error was deemed reasonable given the emphasis in model development on capturing h-current dynamics, including PIR spiking, as discussed in detail in the text.
Indeed, the repetitive spiking frequencies and latencies to the first PIR spike highlighted in Our assertion that this model is appropriate for use in settings beyond those directly constraining   (2006). We note that these differences are also apparent in the experimental data alone, but we 253 focus on the differences in the models given the aims of this study. as shown in Figure 7A (in comparison to experimental results shown in Figure 7B). We note that we 288 will compare these results to those from rodent-derived models in the following section. Given the impetus of this modeling endeavor, we compare the capacity for each of these three 329 models to exhibit subthreshold resonance. In applying an identical ZAP protocol as above for our 330 L5 Human model (see Figure 7), we find that both of these other models, unlike our L5 Human 331 model, exhibit subthreshold resonance at approximately 4.6 Hz as shown in Figure 9. 332 A quantification of these model comparisons is given in Table 2. Alongside results for the 333 baseline models (illustrated in Figure 9), we also include results for the L5 Human and Hay models 334 with all channels besides the h-channel blocked in order to facilitate a more direct comparison with 335 the Kalmbach model (which has no other active ion channels). This alteration results in a minor 336 change in the resting membrane potential (RMP) of the neuron, as would be expected, but no major 337 change in its resonance frequency. 338 The finding that the Hay model exhibits subthreshold resonance is as expected considering 339 that subthreshold resonance has been previously observed in rodent L5 cortical pyramidal cells   and the corresponding peak begin shifting rightwards, with an obvious peak appearing in panels C 356 and D. This resonance is also clearly shown in the corresponding voltage traces. 357 By comparing the different resting voltages in the protocols presented in Figure 10 (and summa-358 rized in Table 3) with the voltage-dependent values in the human h-current model (shown in Figure   359 8), a correlation is apparent between the tendency to exhibit subthreshold resonance and faster   for potentially broader conclusions to be drawn regarding human and rodent differences. 402 Before beginning this investigation, it is important to note that such a switch between human 403 and rodent h-current models would affect other aspects of the cellular model (including, for 404 example, the resting membrane potential, as well as the potential activity of other ion channels) 405 that might affect its behavior. Moreover, the differing morphology and passive properties that 406 make up the "backbones" of these models also differ significantly, and these properties also play a 407 role in dictating a neuron's frequency preference (Hutcheon and Yarom, 2000; Rotstein and Nadim, 408 2014). It is for these reasons that we emphasize that, in performing such a "switch", we create 409 new "hybrid" models that must be approached cautiously. However, a very specific focus on the 410 subthreshold dynamics of these "hybrids" makes their use as presented here reasonable. There 411 are two primary rationales for this assertion: first, a focus on subthreshold dynamics significantly 412 minimizes the role that other ionic currents (whose features vary between "model backbones") will 413 play in the dynamics; and second, by only switching the h-current models (i.e. the kinetics of the 414 h-current), and not the distribution nor conductance of the h-channel, the focus can be mainly on 415 how the different kinetics might play a role (i.e., differences shown in Figure 8). 416 The results obtained are summarized in with -400 pA DC shifts, both the "Hay-L5 Human" and the "Kalmbach-L5 Human" models show model generation we found that a best "fit" to our human experimental data led to significant 546 changes in a variety of conductances (see Table 5) as well as the kinetics of the h-current. Although consisting of a mix of HCN1 and HCN2 subunits display slower kinetics than those seen in HCN1 610 homomeric h-channels (Chen et al., 2001). Taken together, these results provide biophysical support 611 for the hypothesis that the differences in the kinetics of the h-current revealed in this work may be 612 driven by different HCN subunit expression between rodent and human L5 pyramidal neurons. 613 Indeed, in general resonance is a relatively uncommonly observed phenomenon in human 614 neurons (Kalmbach et al., 2018; Chameh et al., 2019), which may be due to the greater expression 615 of HCN2 channels in human neurons generally (Kalmbach et al., 2018). anterior temporal lobectomy for medically-intractable epilepsy (Mansouri et al., 2012). Tissue 647 obtained from surgery was distal to the epileptogenic zone tissue and was thus considered largely 648 unaffected by the neuropathology. We note that this is the same area from which recent data 649 characterizing human L3 cortex was obtained (Kalmbach et al., 2018). 650 Immediately following surgical resection, the cortical block was placed in an ice-cold (approxi- For the subset of experiments designed to assess frequency dependent gain, slices were 667 prepared using the NMDG protective recovery method (Ting et al., 2014)  To supplement the data from the primary neuron, we made use of averaged experimental data 705 from multiple secondary cells. This provided the data of Table 1, which are averaged data from 147 706 cells, and the mean ± standard deviation plots in Figure 6. For the data in Figure 6, we note that for were digitized at 20 kHz using a 1320X digitizer. The access resistance was monitored throughout 727 the recording (typically between 8-25 MΩ), and cells were discarded if access resistance was > 25 728 MΩ. The liquid junction potential was calculated to be 10.8 mV which is corrected for whenever the 729 experimental data is used for modeling or in direct comparison to model values (i.e. Figure 6), but 730 not when the experimental data is presented on its own (i.e. Figure 7B and D).
where is a local constant equal to 2.  Table 5. 838 The locations of each of the 10 ion channel types in our human L5 pyramidal cell model are 840 summarized in Table 6, and utilize a classification of each compartment in the neuron model as part human and rodent cell morphology, even in similar brain regions (Beaulieu-Laroche et al., 2018). 850 As such the region of this increased calcium activity, where the Ca_LVA maximum conductance is 851 multiplied by 100 and the Ca_HVA maximum conductance is multiplied by 10, is chosen to be on 852 the apical dendrite 360 to 600 microns from the soma. 853 The third exception are the h-channels. The function used to model the "exponential distribution" 854 of h-channels along the dendrites (Kole et al., 2006; Ramaswamy where "dist" is the distance from the soma to the midway point of the given compartment, the 859 denominator of "1000" is chosen since this is the approximate distance from the soma to the 860 most distal dendrite, and "gIhbar" is the h-current maximum conductance value that is optimized.

861
"gIhbar" is also the value of the maximum conductance in the soma and basilar dendrites (i.e. the Ih 862 maximum conductance is constant across all compartments in these regions). The first step in the "cycling" model development strategy (schematized in Figure 2)  data obtained from an analogous protocol (Brent, 1976). Here, we fit the model to five different 870 current clamp protocols experimentally obtained from the primary neuron from which we obtained 871 our human L5 cell morphology. As the experimental current clamp data was obtained in the 872 presence of TTX, all sodium conductances were set to zero and not altered in this step. Additionally, 873 the potassium channel currents primarily involved in action potential generation, K_Pst and SKv3_1, 874 were omitted from the optimization and "hand-tuned" in the second step of the cycle. 875 We chose to use three hyperpolarizing current clamp traces, with -400 pA, -350 pA, and -300 pA 876 current amplitudes, because at these hyperpolarized voltages it was reasonable to assume that the 877 h-current was primarily responsible for the voltage changes (Toledo-Rodriguez et al., 2004). This 878 allowed us to accurately fit not only the h-current maximum conductance, but also its kinetics (see 879 Equation 3 above).

880
A hyperpolarizing current step with a small (-50 pA) magnitude was chosen to constrain the 881 passive properties, as near the resting membrane potential it is primarily these properties that 882 dictate the voltage responses ("charging" and "discharging") to a current clamp protocol. We note 883 that this trace does not represent a perfectly "passive" neuron, as some conductances (such as those 884 due to the h-current) are active, albeit minimally, at mildly hyperpolarized voltages (only the sodium 885 channels were directly blocked in this protocol, via the application of TTX). Nonetheless, given 886 that our model fit this current clamp data well, and also mimicked the "charging" and "discharging" 887 portions of all the current clamp protocols included in the optimization, we are confident that we 888 accurately approximated the passive properties of our particular human L5 pyramidal neuron. The 889 final passive properties are shown in Table 5  Finally, a depolarizing current step (100 pA) was chosen to ensure the model was not "overfit" 894 to the hyperpolarized data. Early in the modeling process, we recognized that a "best fit" of the voltage trace in which this channel most affected the voltage, namely the initial "sag" following a 909 hyperpolarizing current steps and the "rebound" in voltage when this inhibition is released. We also 910 chose portions of the voltage trace to emphasize in the error calculation in order to ensure the 911 model cell closely approximated the resting membrane potential observed experimentally, as well 912 as matched the "charging" and "discharging" features heavily influenced by passive properties. 913 We note that these differential "weights" were chosen only after a rigorous exploration of how 914 these choices affected the overall model fit; indeed, this choice yielded a model that both qualita-915 tively and quantitatively best fit the experimentally-observed behavior of our human L5 cortical 916 pyramidal cell. We also note that the possibility that our final parameters represented a "local", 917 rather than "global" minimum in the optimization was investigated by running the optimization with 918 a variety of initial conditions; the solution with the minimum error from all of these trials is the one 919 presented here. 920 Matching of spiking features 921 After optimizing the parameters using MRF, we then tuned the sodium and potassium conductances 922 involved in action potential generation by hand in order to achieve PIR and repetitive spiking 923 behaviors reasonably approximating that seen in experiments (and summarized in Table 1). As 924 described above, we sought to achieve this reasonable behavior while minimizing the relevant 925 potassium conductances so as to best fit the 100 pA current clamp trace. 926 In this step, we also found that a "shift" in the activation curve for Na_Ta (see Equation 4 above) 927 was necessary to achieve PIR spiking as seen experimentally. We sought to minimize this shift for 928 simplicity, but also because a side effect of this leftward shift was an increase in repetitive firing 929 frequency that approached the upper limit of what was biologically reasonable. We note that the The "cycling" mechanism described in detail above was run until there was no significant improve- Moderate constraints were placed on the range of certain parameters in order to ensure that, 959 in finding the best "fit" to the data, these values did not enter a regime known to be biologically 960 unlikely or that would lead to unreasonable spiking characteristics. In order to preserve reasonable 961 spiking behavior, the maximum value for the Ca_LVA maximum conductance was set to 0.001 962 nS/cm 2 , the maximum value for the Ca_HVA maximum conductance was set to 1e-05 nS/cm 2 , and 963 the minimum value of the Im maximum conductance was set to 0.0002 nS/cm 2 . These values were 964 determined after rigorous investigation of the effects of these maximum conductances on the 965 spiking properties. 966 Further constraints were placed on the passive properties of the neuron to make sure the 967 neuron not only matched "charging" and "discharging" properties in the current clamp data, but also The "-20" term in the equation is the "shift" from Kole et al. (2006). The parameters dictating the 1033 model which has non-uniform passive properties and uniformly distributed h-channels (amongst 1034 the soma, apical, and basilar dendrites) are given in Table 7. We ensured our implementation of 1035 this model was appropriate by directly replicating Figure 7B