A population of in silico models identifies the interplay between Nav 1.8 conductance and potassium currents as key in regulating human dorsal root ganglion neuron excitability

Background: The Nav 1.8 sodium channel has a key role in generating repetitive action potentials in nociceptive human dorsal root ganglion neurons. Nav 1.8 is differentiated from other voltage-gated sodium channels by its unusually slow inactivation kinetics and depolarised voltage-dependence of activation. These features are particularly pronounced in the human Nav 1.8 channel and allow the channel to remain active during repolarisation. Gain-of-function mutations in Nav 1.8 have been linked to neuropathic pain and selective blockers of Nav 1.8 have been developed as potential new analgesics. However, it is not well understood how modulating the Nav 1.8 conductance alters neuronal excitability and how this depends on the balance of other ion channels expressed by nociceptive neurons. Methods: To investigate this, we developed a novel computational model of the human dorsal root ganglion neuron and used it to construct a population of models that mimicked inter-neuronal heterogeneity in ionic conductances and action potential morphology Results: By simulating changes to the Nav 1.8 conductance in the population of models, we found that moderately increasing the Nav 1.8 conductance led to increased firing rate, as expected, but increasing Nav 1.8 conductance beyond an inflection point caused firing rate to decrease. We found that the delayed rectifier and M-type potassium conductances were also critical for determining neuronal excitability. In particular, altering the delayed rectifier potassium conductance shifted the position of the Nav 1.8 inflection point and therefore the relationship between Nav 1.8 conductance and firing rate. Conclusions: Our results suggest that the effects of modulating Nav 1.8 in a nociceptive neuron can depend significantly on other conductances, particularly potassium conductances.


Introduction
The sodium channel Nav 1.8 is selectively expressed in small diameter dorsal root ganglion (DRG) neurons and has an important and unique role in pain signalling, due to its slow and depolarised kinetics of inactivation [Akopian et al., 1996;Sangameswaran et al., 1996] and fast recovery from inactivation [Dib-Hajj et al., 1997] compared with other voltagegated sodium channels.Human Nav 1.8 has particularly slow inactivation and larger persistent current compared to rodent Nav 1.8 [Han et al., 2015].Gain-of-function mutations in Nav 1.8 contribute to painful peripheral neuropathy [Faber et al., 2012;Huang et al., 2013;Han et al., 2014], and Nav 1.8 has been shown to support repetitive firing and the characteristic broad shoulder of DRG neuron action potentials [Renganathan et al., 2001;Blair and Bean, 2002].Nav 1.8 has also been identified as a potential therapeutic target, due to its effects on DRG neuron excitability, and selective blockers of the channel have been developed [Jarvis et al., 2007;Payne et al., 2015].However, to date, there has only been one clinical trial of a selective Nav 1.8 blocker, which was an unsuccessful trial of PF-04531083 [Payne et al., 2015] for post-surgical dental pain.
It is therefore important to better understand how Nav 1.8 affects the excitability of human DRG (hDRG) neurons.However, the effect of an ion channel on neuronal electrophysiology depends not only on its individual kinetics and expression level but also on the background ensemble of other ion channels expressed in the neuron.For example, in DRG neurons expressing Nav 1.8, a gain-of-function Nav 1.7 mutation caused hyper-excitability while in sympathetic ganglion neurons without Nav 1.8 the same mutation causes hypo-excitability [Rush et al., 2006].Understanding how variation in multiple ionic currents will interact with each other, not just on a qualitative level (whether a channel is expressed or not) but on a quantitative level (the level of expression of each channel), is challenging.This is an important consideration though, as ion channel densities have been shown to have high inter-neuronal variability [Schulz et al., 2006].They can cause different neurons to respond differently to the same experimental conditions, e.g.application of the same concentration of a drug.Pathological conditions also substantially alter expression of multiple channels in DRG neurons [Black et al., 2004;Cao et al., 2010;Zhao et al., 2013].
In this study, we investigate the effects of modulating the Nav 1.8 conductance (GNav 1.8) on hDRG neuron excitability using computational modelling and simulation.We hypothesised that altering the Nav 1.8 conductance would be positively correlated with firing rate, but that the quantitative change in firing rate would depend on the balance of other conductances in a particular model.To test this hypothesis, we develop a novel model of the hDRG neuron based on human-specific ionic currents.We construct a population of 328 experimentally calibrated human DRG neuron models [Britton et al., 2013;Marder and Taylor, 2011;Prinz et al., 2003/4] that integrate data from electrophysiological recordings of human DRG neurons at the cellular [Davidson et al., 2014] and ion channel levels [Han et al., 2015;Vasylyev et al., 2014;Huang et al., 2014].Each model in the population shares the same underlying kinetics and equations for the eight ionic currents in the model, but has a different set of ionic conductances, reflecting inter-neuronal heterogeneity in ion channel densities.Importantly, every model in the population produces action potential behaviour in range with the experimental variability reported in human DRG neurons by Davidson et al. [2014] when simulated with equivalent protocols and measured with eight action potential parameters.
This study advances the replacement of animal DRG neurons for research into human DRG neuron electrophysiology by developing the first computational model of human DRG neuron electrophysiology and providing open-source software to run simulations using the model.This software allows interested users to begin using the model with the population of models approach described in this study for modelling inter-neuronal variability.The main situations where we think our model would be particularly suited to replacing the animal DRG neuron experiments are for hypothesis generation, as many different experimental conditions can be rapidly simulated and analysed; and for detailed analysis of ionic current behaviours as each current can be output and visualised.

Overview of modelling approach
In this study, we use a non-standard modelling approach based on experimentally calibrated populations of models [Britton et al., 2013;Marder and Taylor, 2011].This methodology, which originated in neuroscience [Prinz et al., 2003[Prinz et al., , 2004;;Marder and Taylor, 2011], has also been used extensively to model the electrophysiology of cardiac cells [Britton et al., 2013;Passini et al., 2017], but to our knowledge this is the first time it has been used to study DRG neurons.This approach allows us to integrate the observed experimental variability in hDRG neuron action potential (AP) parameters [Davidson et al. 2014] into the modelling process, instead of the standard approach where a single model parameterised to reflect a typical neuron is used [Marder and Taylor, 2011].
Here we explain the population of models methodology.For more detail, see [Britton et al., 2013;Marder and Taylor, 2011;Muszkiewicz et al., 2016].The main idea is to create a group of neuronal models that all have different underlying conductance values (representing inter-neuronal variability in ion channel densities) that, when simulated, all generate voltage traces with AP biomarkers that are within the range of values observed experimentally in hDRG neurons.
To create our population of models, we first construct a baseline model.This model is a Hodgkin-Huxley-type neuronal AP model [Hodgkin and Huxley, 1952] comprising a set of equations describing the changes in membrane voltage and ionic current gating variables over time.We then select parameters in the model to be varied to represent inter-neuronal variability.In this study we make the assumption that ion channel expression varies from neuron to neuron but that the structure of each ion channel type does not.Therefore, we randomly vary all conductances in the baseline model but leave the kinetics parameters of each ionic current constant.Through random variation of all conductances we obtain a pool of candidate models, each with a different set of conductances.We simulate each model to find its rheobase, record its voltage trace and from this trace calculate AP biomarkers.We then calibrate the population by comparing each AP biomarker in each model to the experimental range for that parameter.If any AP biomarker for a model is outside the experimental range, the model is rejected.If all tested AP biomarkers are within their experimental ranges then the model is accepted.The final population of models consists of only the accepted models.All models in the population have a different set of conductances but each model is in range with experimental data on AP parameters of ex vivo hDRG neurons.

Baseline model
For this study, we chose to model the human dorsal root ganglion neuron soma, rather than the axon or full DRG neuron geometry.This is because the soma was the section of the neuron that was patched and recorded from in both the experimental data we used to calibrate the population [Davidson et al., 2014] and in the vast majority of DRG neuron electrophysiological studies, due to its large size compared with the axon.Our aim was to create a model that could be directly compared against results from the experimental system from which data are available.
The kinetics (voltage dependencies of activation and inactivation, and time constants) of each ionic current were not varied in this study, while all eight conductances were simultaneously varied.These choices are based on the assumptions that: 1) kinetic parameters depend primarily on ion channel structure, which as a simplifying assumption for modelling we assume will not vary between different DRG neurons; 2) conductances represent the density of each ion channel type in the cell membrane, and these vary substantially between neurons [Schulz et al., 2006].
The soma was modelled as a single cylindrical compartment with length=30 μm, diameter=46 μm, and cytoplasmic resistivity=1 Ω/cm following the values used in the DRG modelling study by Choi and Waxman [Choi and Waxman, 2011].A cylindrical geometry was used in line with how the NEURON modelling software represents all neuronal compartments.The specific membrane capacitance was set to the commonly used value of 1 μF/cm 2 .The internal and external ionic concentrations were fixed to match the experimental conditions used by Davidson et al. [2014] with extracellular [Na + ] =145 mM, intracellular [Na + ]=5 mM, extracellular [K + ]=3 mM, and intracellular [K + ]=135 mM.
The equations for the baseline model can be found in the supplementary material and source code is provided in the data supplement.

Experimental data and stimulus parameters
We used data on experimental AP biomarkers of ex vivo human DRG neurons from Davidson et al. [2014] to calibrate the population of models against experimental data.This dataset included data from recordings of 141 ex vivo human DRG neurons from five organ donors.Each neuron's soma was recorded under two stimulus protocols, an 800-ms step stimulus and a 500-ms ramp stimulus.In both cases, stimulus current was increased in steps of 50 to 100 pA until AP threshold was reached for that neuron, and AP biomarkers were calculated from traces at that stimulus current amplitude.

Measuring AP properties
We used eight AP biomarkers to quantify the electrophysiological properties of hDRG neuron APs: AP threshold voltage; AP peak voltage; AP slope maximum; AP slope minimum; AP full width; after-hyperpolarisation time constant; resting membrane potential; and rheobase.Each biomarker was calculated as follows:

Rheobase
Rheobase was calculated for each model by running step current simulations with increasing stimulus current amplitudes in increments of 100 pA from 0 pA to up to 5000 pA until the first simulation in which one or more action potentials were detected.The stimulus amplitude of that simulation is the rheobase for that model.

AP threshold voltage
Following the approach used by Davidson et al. [2014], AP threshold voltage was calculated from a ramp current simulation and defined as the voltage during the upstroke of an action potential at which the voltage-time gradient first surpasses 5 mV/ms.

AP peak voltage
AP peak voltage was calculated as the peak voltage attained during an action potential.

AP slope maximum and minimum
AP slope minimum and maximum was calculated as the maximum voltage-time gradient and minimum voltage-time gradient during an action potential.

AP full width
AP full width for an action potential was calculated as the duration between the time at which the voltage first crossed the AP threshold voltage from below, and the time at which the voltage first crossed the AP threshold voltage from above.

After-hyperpolarisation time constant
The time constant (τ) of the after-hyperpolarisation of an action potential was calculated by fitting the afterhyperpolarisation to the following single-exponential model with a curve fitting algorithm and extracting the value of τ: The after-hyperpolarisation of an action potential was defined as the period from the minimum voltage that occurred after the AP peak, up to the point at which the voltage gradient exceeded 5 mV/ms or 50 ms had elapsed, whichever came first.

Resting membrane potential
The resting membrane potential of a model was defined as the minimum voltage that occurred during the last 10% of a simulation in which no stimulus current was applied.

Sampling and calibrating the population of models
To create the initial pool of candidate models, all conductances in the baseline model were simultaneously and randomly sampled in a range of 0-2-times the baseline model value for that conductance (Extended data -Equations_and_ conductance_densities.docx-Conductance densities table) to create 20,000 candidate models with the same ionic currents and equations but different conductance values.Conductance sampling was performed using Latin hypercube sampling (McKay et al. 1979), which is a technique that evenly samples each conductance without exhaustively sampling all possible conductance combinations.We used Latin hypercube sampling as sampling all possible combinations of eight conductances and simulating the resulting models is not computationally feasible, even if each conductance is sampled coarsely.The range of 0-2-times baseline was chosen to allow a broad range of different conductance profiles within the population.Previous studies of neurons have shown high variability in mRNA expression and current densities between healthy neurons of the same type [Schulz et al., 2006].We set the minimum of our sampling range to zero (equivalent to no expression of functional ion channels) to allow our population to potentially include model neurons that produce a normal AP under normal conditions even without substantial expression of all ionic currents in the model.This can represent conditions such as a loss of function mutation that may not affect normal neuronal function but could change the response when conditions are changed, e.g. if another channel is blocked by a drug.
Each model was simulated under 800-ms step and 500-ms ramp protocols at its individual rheobase to match the stimulation protocols used in Davidson et al. [2014].From the voltage traces from these simulations we extracted values for the eight AP biomarkers described previously: AP threshold voltage; AP peak voltage; AP slope maximum; AP slope minimum; AP full width; after-hyperpolarisation time constant; resting membrane potential; and rheobase.
We then calibrated our candidate models against experimental ranges of these AP biomarkers calculated from data in Table 3 of Davidson et al. [2014].The experimental range for each AP biomarker was defined as the range spanned by the mean AE1.5 standard deviations of the data for that parameter.The choice of 1.5 standard deviations was partly arbitrary but was chosen so that the range would include approximately 90% of the total probability density of the distribution of each AP biomarker assuming a normal distribution.This was so the population would capture a large proportion of total variability in each AP biomarker while excluding extreme outliers.

Varying conductance parameters in simulation with scaling factors
In this study, we performed simulations to investigate the effects of altering the GNav 1.8 conductance, as well as the GKdr and GKM conductances, across the whole population.In each of these simulations, we implemented this by multiplying the conductance of each model in the population by a fixed scaling factor.This mimics how drug block would affect the conductances of different neurons.For example, two neurons with GNav 1.8 values of 1 and 2 respectively subject to a scaling factor of 0.5 would have their GNav 1.8 values reduced to 0.5 and 1, respectively, equivalent to 50% Nav 1.8 channel block.

Simulations and software
All simulations were run using NEURON version 7.5 (RRID:SCR_005393) [Carnevale, 2006], through NEURON's Python (RRID:SCR_008394) interface.NEURON uses the adaptive time step ODE solver CVODE [Hindmarsh et al., 2005] to solve model equations.All simulations were recorded at 40 kHZ (NEURON's default) then down-sampled to 20 kHz to match the sampling rate used by Davidson et al. [2014].All analysis and visualisation was performed using Python 3.7.The Python module 'drgpom' we developed to run NEURON simulations on populations of models and analyse the results is included in the extended data and the most up to date version can also be found at [https://github.com/oliverbritton/drg-pom].The module is also available on the Python online package repository Pypi [https://pypi.org/] and so can be downloaded and installed with a single command.The module comes with a series of examples to guide users in how to use the module, which we hope will help potential users, who have not used populations of models before, run their own simulations.
Visualization of balances of ionic currents during simulations using currentscapes An advantage of computer simulations is the ability to simultaneously record the voltage trace and all individual ionic currents in each simulation.To visualise the relative contributions of all ionic currents throughout simulations we used a visualization technique called a currentscape [Alonso and Marder, 2019].A currentscape shows the voltage trace and the fractional contribution of each ionic current throughout a current clamp simulation.Each currentscape shows the voltage trace from a simulation at the top, followed by plots showing the fractions of each outward current in the model and each inward current in the model.This allows us to see which currents are influential at each phase of the action potential, and how this differs from model to model.

Protocols for using and installing drgpom
This guide for installing the drgpom package for Python assumes you are installing on a Windows machine and has been tested on Windows 7 and Windows 10, and should also work on Linux (drgpom has been tested on the Fedora Linux distribution).

Prerequisites
The use of this software requires a basic familiarity with Python.One of the best places for scientists to learn Python for scientific computing is to follow the lessons from the Software Carpentry website: https://software-carpentry.org/.It also requires installation of Python 3 and the NEURON simulation software, both of which are freely available and open source.
Installing the dependencies for drgpom 1. Download and install the Anaconda Python 3.7 distribution (https://www.anaconda.com/).This is a Python distribution for scientific computing including most of the key libraries needed for drgpom.
2. Install NEURON (https://neuron.yale.edu/neuron/).NEURON is simulation software that is widely used in computational neuroscience for simulating neuronal electrophysiology.We use it to run the individual simulations and implement the models in drgpom.

Installing drgpom and compiling models
The simplest way to install drgpom is to download it using pip, which is a Python tool for automatically download and installing packages from online repositories.If you have installed Anaconda command window, pip will come installed with it.
You can install drgpom using pip by opening a command window (in Windows this will be a command prompt and can be opened by running the command cmd from the Run option in the Start Menu, in Linux this will be a terminal).From there run the command: This should produce some output including the message "This simulation has finished running" and the results of the simulation along with simulation metadata such as model details and parameter values will be saved in a file called "example_population.pkl".This is a pickle file, which is a format for saving Python objects so they can be reloaded.A set of images will also be saved showing the voltage traces from each model in each simulation.To see how to view results from this simulation see the example notebook.
Configuring your own simulations using a pre-existing population of models 1.To run a new simulation using the provided population of models, copy "example_simulation.py" and alter the options parameters to configure your simulation.The simulation parameters for this type of simulation are described in Table 1 below.2. To run your simulation, open a command window or Jupyter notebook in the directory the simulation script is saved in.
3. If using a command window type: python script_name.pycores=n to run the script, where n is the number of CPU cores you want to use (or don't include cores=n to use all but one core by default).If using Jupyter notebook, run the code %run script_name.pycores=n.In both cases the simulation will run, leave the notebook or prompt running until it has finished.When the simulation has finished you will see a message that looks like:

À65
flags Special flags for particular hardcoded options An empty dictionary for a dictionary or flag name and flag value pairs (see Simulation 2 in example_simulation.pyfor an example).

{} (empty dictionary)
This parameter is mainly used to allow AP full width for step stimulus simulations to be computed using the threshold voltage from a corresponding ramp stimulus simulation.where your_simulation.py is the name of your simulation file.
3. The simulation will now run, do not close the window before it has finished as this will halt the simulation.
4. Once the simulation is completed the results will be saved to a file that can be opened in Python.The examples notebook provides example code showing how to view your results.

Notes on the model_details parameter
The model_details parameter defines the ion channel models that will be included within each model in the population, and which parameter(s) in each ion channel model will be varied.To define this parameter we first initialise a dictionary called model_details containing another dictionary called mechanisms that will store the details of each ion channel model: 'param_name' is the name drg-pom will know that parameter by, for example 'GNav 1.7'.You can choose this parameter name to be whatever you want, ideally something clear and memorable.These parameter names should correspond to those used in parameter_data to define the scaling ranges each parameter will be sampled over.
'param_name_in_nrn' is the name that parameter is known by in NEURON.The format of this name is the name of that parameter in the.mod file of its associated ionic current model, followed by an underscore and nrn_mech_name.For example conductances are generally given the variable name 'gbar', so the conductance for the kleak ion channel model is named 'gbar_kleak'.See the example_population_creation.pyfile for an example of a complete model_details parameter.

Running a simulation with multiple channel block conditions on multiple channels
The simulation script vary_amp_and_gnav18.py in the examples directory provides an example of a simulation design similar to those used in the paper where multiple simulations are run simulating changes in multiple conductances at multiple different levels.This example can be used as a template to design your own simulation studies.

Properties of the population of human DRG neuron AP models
The population of hDRG neuron AP models consisted of 328 models each of which had the same underlying equations but a different combination of eight ionic conductances.All models in the population had substantially different combinations of conductances, but all of them produced AP biomarkers at rheobase that were within 1.5 standard deviations of the mean values of each parameter, based on data recorded from ex vivo hDRG neurons [Davidson et al., 2014].
Voltages traces from the population of models (Figure 1A) show they reproduce key characteristics of the hDRG neuron AP, including its characteristic broad shoulder [Davidson et al. 2014, Han et al., 2015], the slow rise to threshold followed by a rapid upstroke, and the presence of an after-hyperpolarization following repolarization.A schematic of the baseline model is shown in Figure 1B.
The traces in Figure 1A were produced by models with a diverse range of conductance profiles (Figure 2).Most conductances (except for GNav 1.8 and GKA) span the sampled conductance range.For GNav 1.8 and GKA, the accepted models only have a subset of the sampled conductance range (GNav 1.8: 0.15 -1.04, GKA: 0 -1.22).For GNav 1.8, this is because models with low GNav 1.8 have peak voltages that are below the experimental range, and models with high GNav 1.8 have threshold voltages that are below the experimental range.For GKA, no model with high GKA had a combination of threshold voltage, peak voltage, and after-hyperpolarisation time constant that was simultaneously within the experimental ranges of these three AP biomarkers.
Most of the pairs of conductances are uncorrelated (Figure 2), however GNav 1.8 and GKdr have a strong positive correlation (r = 0.82).This correlation is caused by calibrating on a combination of three biomarkers: threshold voltage, peak voltage and AP full width parameters, as removing all of these AP biomarkers from the calibration process also removes the correlation between GNav 1.8 and GKdr and removing any one of them reduces the correlation coefficient.This indicates that in our model the balance of these two currents is important throughout the action potential, from upstroke (threshold, peak voltage) to repolarisation (AP full width).

Visualising the balance of currents in different hDRG models
As can be seen in Figure 2, we found many different combinations of conductances that produced viable models consistent with observed experimental variability in hDRG AP biomarkers.However, while each of these models produced similar voltage traces, the different conductances in each model mean that the balance of underlying ionic currents in each model were substantially different.Figure 3 shows the conductances and relative magnitudes of each current in six models from the population using a visualisation technique called a currentscape [Alonso and Marder, 2019] (see Methods for details).
We see that during quiescence the primary active currents are IKdr (purple) and IKleak (pink) for the outward currents and Ih (brown) and INav 1.9 (green) for the inward currents.The relative proportions of these currents are however highly variable from model to model.Similarly, the main currents during the upstroke are INav 1.7 (red) and INav 1.8 (blue), but some models have a large contribution by Nav 1.7 and in some models the contribution is very small, while the contribution of INav 1.8 is uniformly large.
During the body of the AP, in all models IKdr (purple) is the dominant outward current and INav 1.8 (blue) is the dominant inward current.These two currents appear to have unique roles in shaping the AP that cannot be performed by any of the other currents in our model.This can explain why GKdr and GNav 1.8 are highly correlated in Figure 2; the values of both conductances must be balanced against one another as no other currents can perform their roles in shaping the AP.Following repolarisation, the potassium currents IKA (orange) and IKM (yellow) are both important for shaping the afterhyperpolarization.However, the balance of these two currents differs from model to model.For example, the model in Figure 3E has most of its post-repolarisation outward current conducted through IKA with a very small contribution from IKM, while the model in Figure 3F relies primarily on IKM with a small contribution by IKA.As we show later, the balance of these currents is not as critical as the balance of GNav 1.8 and GKdr at rheobase, where only a single AP is fired.However, at higher stimulus amplitudes that support repetitive firing, the value of GKM becomes important for modulating firing rate.
The Nav 1.8 conductance has non-linear effects on firing rate Following the evaluation of the population of hDRG neuron models shown in the previous section, we investigated how changing the Nav 1.8 conductance affects hDRG neuron excitability across a wide range of different backgrounds of other ionic conductances.We therefore simulated different levels of block and enhancement of the Nav 1.8 conductance in every model in the population, across a range of stimulus amplitudes. Figure 4A shows the GNav 1.8 scaling factor and stimulus amplitude parameters for the 441 simulations that we ran.For each simulation the 328 models in the population  were each simulated under an 800 ms step stimulus current protocol.The Nav 1.8 conductance was varied by multiplying GNav 1.8 in each model by a scaling factor, which is similar to how drug block or enhancement would affect the conductance.
Figure 4B shows the mean firing rate (APs/s) of the population of models in each simulation.Mean firing rate is not maximised at the highest GNav 1.8 we simulated but is instead maximised when GNav 1.8 is slightly above the baseline value (between 1.1 and 1.3 times baseline, depending on the stimulus amplitude of the simulation).If GNav 1.8 is increased beyond this point, the mean firing rate decreases.This result implies that increased expression of Nav 1.8 could lead to a decrease in firing rate.We also find that reducing GNav 1.8 to 0.5 times baseline or less is sufficient to stop repetitive firing at all of the stimulus amplitudes we simulated, indicating the importance of Nav 1.8 in our model for enabling repetitive firing.
We hypothesised that one reason for the decrease in firing rate as GNav 1.8 is increased could be increased AP width due to increased Nav 1.8 current during AP repolarisation.A wider AP could reduce the maximum firing rate of a neuron by lengthening the minimum time between subsequent APs.Unlike other sodium channels Nav 1.8 is active during AP repolarisation, contributes to the characteristic shoulder of the DRG action potential [Han et al., 2015], and is therefore important for determining AP duration.Figure 4C shows that mean AP half width (defined as the period during an AP in which voltage is continuously above a threshold halfway between the threshold voltage and peak voltage) increases steadily with increased GNav 1.8 and is insensitive to stimulus amplitude.For example, averaging over different stimulus amplitudes, mean half-width of the population changed from 3.6 ms at baseline (GNav 1.8 = 1) to 4.5 ms when GNav 1.8 is doubled (GNav 1.8 = 2).
Plotting the relationship between half width and firing rate in individual models (Figure 5) shows that firing rates within the physiological range (approximately 0 -25 AP/s) can occur across a range of half-widths, but the very high firing rates that some models in the population exhibit (up to approximately 120 AP/s) only occur at the lower half-widths found predominantly in simulations with lower GNav 1.8 scaling values (e.g. Figure 5, second panel from the top).However, at the lowest values of GNav 1.8, there is very little repetitive firing (Figure 5, top panel).Therefore, one of the effects of Nav 1.8 on hDRG neuron electrophysiology may be to support repetitive firing but also limit the maximum firing rate by broadening the AP.
Very high firing rate models are distinguished from other models by low GKM and high GKDr potassium conductance We noticed that a substantial fraction of models fired at a very high rate (> 20 AP/s) (Figure 5.This is above the physiological values reported from real hDRG neurons [Meyer and Campbell, 1981] but these models are still relevant because they exhibit a normal hDRG-like AP at rheobase.The ways that these models' conductance profiles differ from models in the population that do not fire at a very high rate could indicate which currents act as brakes or enhancers of excitability when the Nav 1.8 conductance is altered. Therefore we divided the population of models into two groups based on the number of simulations from Figure 4A in which a model fired at a rate > 20 AP/s.The top 25% of models were assigned to the "most often rapidly firing" group" and the other 75% of models to the "all other models" group.Figure 6A shows the distribution of each conductance in each group.Interestingly, GKdr and GKM differ substantially between the two groups.GKM is much smaller in the most often rapidly firing group (0.64 AE 0.32) than the other models group (1.41 AE 0.43).GKdr is larger in the rapidly firing group (1.60 AE 0.44) compared to the other models group (1.28 AE 0.43).That GKdr is actually larger in the rapidly firing group suggests it is not just a lack of total potassium current conductance, and a corresponding lack of outward current to oppose repetitive firing, that causes these models to fire so rapidly.
We therefore simulated the effects of modulating each combination of Nav 1.8 and one of the other seven conductances in the model to confirm whether the firing rate of models in the population were sensitive to GKdr and GKM, but no other conductances.Modulating either GKdr or GKM substantially altered the effects of varying GNav 1.8 on firing rate (Figure 6B), while the other five conductances showed little effect on firing rate (Figure 7).Modulating GKdr (Figure 6B, top row) altered firing frequency by increasing the Nav 1.8 values required for high frequency repetitive firing as GKdr was increased.Modulating GKM primarily reduced firing rate across all GNav 1.8 values as GKM was increased.
M-type and delayed-rectifier potassium currents modulate the effects of Nav 1.8 and collectively control excitability As we found that GKM and GKdr levels each have an important role in determining the excitability of models in the population (Figure 6B), we ran a simulation study where we co-varied four parameters: stimulus amplitude; GNav 1.8; GKdr; and GKM, to see how changes in these three important conductances interacted.As we were simulating changes in four parameters, we used a coarser range of values for each parameter compared to the previous simulation study design (Figure 4A).Stimulus amplitude was varied from 0 to 6000 pA in increments of 1000 pA and the conductance scaling factors for GNav 1.8, GKdr and GKM were each varied from 0.0 to 2.0 in increments of 0.5.All combinations of the four parameters over these ranges were simulated for a total of 875 simulations on each model in the population.
Figure 8 shows the mean firing rate of the population in each of these simulations.As plots go from left to right, GKM increases, as plots go from bottom to top, GKdr increases.To better understand the relationship between the values of GKdr and GKM and their effects on firing rate, we looked at the overall trends for these conductances.We averaged mean firing rate over all of the GNav 1.8 scaling values we simulated.This condenses the information in each panel of Figure 8 down to a single averaged firing rate.We then plotted curves of either GKM or GKdr against mean firing rate averaged over GNav 1.8.This produced a group of curves (Figure 9) that show the general trend from varying GKdr (Figure 9A) or GKM (Figure 9B).
Broadly we see three different relationships between firing rate and the scaling factor of GKdr and GKM, depending on the scaling factor of the other conductance: 1. Monotonically negative response -increasing the potassium conductance that wasn't fixed always decreased firing ratethe expected response for a potassium current.This occurred when GKM or GKdr scaling factors were fixed at 1.5 and 2 (Figure 9A and B, red and purple curves).
2. Monotonically positive response -increasing the potassium conductance that wasn't fixed always increased firing rate.This is opposite to the expected response for a potassium conductance and occurs because without sufficient potassium current the models cannot effectively repolarize, which prevents repetitive firing.This only occurred when GKM was fixed at 0.
3. Inflected responseincreasing the potassium conductance that wasn't fixed increased firing rate up to an inflection point, beyond which increasing the conductance further caused firing rate to decrease.This is the most common scenario.
These results shows that in our model, across a wide range of ionic conductance profiles, the response of firing rate to a change in a potassium conductance can qualitatively differ depending on the value of another potassium conductance.Therefore channel block by a drug could increase or decrease firing rate, depending on the ionic background of the hDRG neuron the drug was applied to.

Main findings
In this study, we investigated the roles the voltage-gated sodium channel Nav 1.8 and the IKdr and IKM potassium currents play in determining hDRG neuron excitability using a population of computational models calibrated against experimentally determined hDRG AP biomarker ranges.We used this population of models methodology to allow us to simulate the variability in ion channel densities found in real neurons and investigate how variability in many different ionic currents collectively interact to determine overall cellular excitability.
The main findings of this study were: 1. We developed a new model of the hDRG neuron AP using published data from hDRG neurons at the AP level and at the ion channel level where data was available.We used this model as a baseline to develop a population of models that integrated experimentally observed variability in hDRG AP biomarkers.
2. Altering the Nav 1.8 conductance had nonlinear effects on excitability.Decreasing the Nav 1.8 conductance from baseline in each model decreased the mean firing rate of the population, however increasing the Nav 1.8 conductance initially increased firing rate up to a point, beyond which firing rate decreased.
3. Simulation studies with the population of models identified key roles for the M-type potassium current and delayed-rectifier potassium current.Increasing the M-type potassium current conductance decreased firing rate.However, increasing the delayed-rectifier potassium current increased firing rate, but only when Nav 1.8 was also increased by a similar factor.
4. Several of our results agree with the literature on DRG neuron electrophysiology, which provides credibility for our modelling approach.Increasing Nav 1.8 conductance increased AP width, in line with current understanding of the role of human Nav 1.8.[Han et al., 2015].Our finding that increasing M-type potassium channel conductance decreased excitability agrees with studies that show M-channel opening drugs such as retigabine reduce sensitivity to nociceptive stimuli in peripheral nerve fibres [Passmore et al. 2012, Vetter et al. 2013, Huang et al. 2016].

The hDRG neuron model
The baseline hDRG model we developed has several novel features that make it adapted to modelling hDRG neurons.We used ionic current models developed from recordings of human channels, for Nav 1.7, Nav 1.8 and Nav 1.9 [Vasylyev et al., 2014;Han et al., 2015;Huang et al., 2014].It was particularly important to have a human-specific Nav 1.8 model as it has substantial differences in kinetics compared to rat Nav 1.8 [Han et al., 2015].The kinetics of inactivation of human Nav 1.8 are slower, and the V1/2 of inactivation is more depolarized, so more human Nav 1.8 channels are available, for longer, during the action potential, compared to rat Nav 1.8.Previous models and most studies of DRG neuron electrophysiology have used rodent DRG neurons due to the difficulty of acquiring hDRG neurons for research.Therefore, species differences in AP morphology and ion channel behaviours [Han et al. 2015] were not captured.
To mimic inter-neuronal variability in ion channel conductances and AP biomarkers we used the population of models methodology [Britton et al., 2013] and experimental data from ex vivo hDRG AP recordings to calibrate the population of models to match the ranges of AP biomarkers from those experiments.Similar approaches to this have been used in many modelling studies of other neuron types (e.g.[Prinz et al., 2003[Prinz et al., , 2004]]) and of cardiac cells (e.g.[Britton et al. 2013, Passini et al., 2017, Muszkiewicz et al., 2016]) but to our knowledge they have not previously been used in studies of DRG neurons.
Where hDRG neuron-specific data was not available we built on previous modelling studies of DRG neurons, in particular work by Tigerholm et al. (2015) and Choi and Waxman (2011).Tigerholm et al. investigated activitydependent slowing in a model of a DRG neuron axon.As human-specific potassium current and hyperpolarisationactivated cation current models were not available we used the same potassium and h-current models that were used in this study.Choi and Waxman investigated how varying GNav 1.7 and GNav 1.8 affected excitability in a single compartment model of a DRG neuron cell body.We adapted the geometry and passive membrane properties used in this study for our model.

Implications for the 3Rs (reduction, refinement and replacement) of animals in research
Computational hDRG neuron models, along with other alternatives such as ex-vivo [Davidson et al., 2014] and stem cellderived hDRG neurons, offer alternatives that could replace the use of rodent DRG neuron recordings in the future.For example, a computational model could be used to investigate how different ion channel modulating drugs interacted with a pathological ion channel mutation in hDRG neurons, to determine if they could reduce or counteract the pathological effects of the mutation.This type of study would require a large number of rodent DRG neurons if performed using an animal model, as the experiment would need to be repeated for each different drug tested.Computational models also have the advantage that they can be refined and iterated upon as new data becomes available, to make the most of the limited availability of hDRG neurons and recordings.Furthermore the fact that they are human-based allows for easier clinical translation, as animal moedels present patho-physiological differences with human cells, which can impact therapy outcomes.
To estimate the number of animals currently used in studies of rodent DRG neuron excitability, that could potentially be replaced by these alternatives, we searched PubMed to estimate the use of rodents in this class of experiment.Our base search query was: (mouse OR mice OR rat OR rodent) AND (nociceptor OR DRG OR dorsal root ganglion) AND pain AND (ion channel OR ion channels OR Nav) with the addition of: 1) AND (drug block OR channel block OR pharmacological block OR current inhibition OR current block) 2) AND (mutation OR channelopathy) These queries aimed to find papers that investigated: 1) Drug-induced changes in human DRG excitability and 2) Changes in human DRG excitability caused by ion channel mutations, as these are the main areas that could be investigated with our population of models.
The base search produced 1018 papers published over the last 5 years (2015-20), indicating wide use of these models.
Filtering further using Queries 1 and 2 reduced these papers to 361 unique non-review papers published between 2015 and 2020.We analysed the first 20 papers for each query to determine average numbers of animals used for DRG studies only.Mean animals used for DRG studies in the 16/40 papers where number of animals used was reported or could be inferred from figures was 103, with a range of 24 to 273.Therefore, annual animal use in these experiments is estimated at 7,400 animals (103 animals multiplier by 361 papers divided by 5 years), while total annual animal use in all DRG experiments is estimated at 21,000 animals (using the figure of 1018 papers found by the base query).These estimates do not include animal use in industry, which is likely to be substantial given the current interest in developing selective sodium channel blockers targeting human DRGs.

Nav 1.8 and GKdr
Changes to ionic currents rarely affect only a single channel.For example, pathological changes in channel expression often affect multiple channels [Black et al., 2004;Cao et al., 2010;Zhao et al., 2013], and even selective ion channel blocking drugs often have non-negligible side effects.
Nav 1.8 is known to support repetitive firing in DRG neurons [Renganathan et al., 2001;Blair and Bean, 2002], due to its rapid recovery from inactivation [Dib-Hajj et al., 1997].Therefore we expected that increasing GNav 1.8 would likely increase firing rate across the population.However, our simulations showed a more complex relationship.Increasing GNav 1.8 by up to approximately 30% increased mean population firing rate but beyond this inflection point increasing GNav 1.8 caused the mean population firing rate to decrease.We also found that the GNav 1.8 value at which this inflection point occurred was highly dependent on GKdr (Figure 9B).Enhancing or blocking GKdr caused the inflection point to move to higher or lower GNav 1.8 values respectively.
The opposing contributions of Nav 1.8 and IKdr during repolarisation of the action potential (Figure 3) suggest the balance between Nav 1.8 and IKdr may be important for supporting repetitive firing in hDRG neurons.However, as more ionic currents are expressed in real hDRG neurons than are included in our model, it is likely the true relationship will be more complicated than simply balancing GKdr against GNav 1.8.

M-type potassium current
We found that the conductance of the M-type potassium current, primarily carried by Kv 7.2, 7.3 and 7.5 in DRG neurons [Passmore et al., 2003], strongly modulated the excitability of the population of models.Increased GKM decreased mean firing rate of the population of models across a wide range of stimulus amplitudes and GNav 1.8 and GKdr scaling factors (Figure 8).This is in agreement with studies showing that M-current activators such as retigabine and flupirtine can reduce sensitivity of nociceptors to nociceptive stimuli [Passmore et al., 2012;Vetter et al., 2013;Huang et al., 2016].The accepted role of the M-current is to stabilize the membrane potential against small depolarisations.Our results suggest that the M-current could also have a role in opposing large stimulus currents during sustained firing, due to the M-current's slow kinetics of activation that lead to it reaching its maximum amplitude only after a period of repetitive firing (Figure 3).

Future development and limitations
We think that future development of the model presented here should focus on three areas: validation against new experimental data as they become available; incorporating additional ionic currents, particularly calcium-activated potassium currents and calcium currents; and expanding the geometry of the model to include the DRG neuron axon and t-junction with the cell body.
The population of models used in this study was calibrated using AP biomarker data from wild type hDRG neurons stimulated under standard conditions (e.g.no applied drugs).Adding additional calibration steps with AP biomarkers from experiments under very different conditions (e.g.under drug block) would further constrain the allowed conductance profiles of accepted models which should result in a population that behaves more realistically under a wider range of simulated conditions.This overcomes an important limitation from previous models, which are based in animal data.
Our current baseline model contains eight ionic current models that cover a diverse range of sodium and potassium currents but do not include calcium currents, calcium-activated potassium currents, or ionic pumps and exchangers.All of these would be good targets to add to the model however there is less data, particularly DRG-specific data, to parameterise models of these currents, so this work would either require using models developed for other neuron types or new experimental data.
Extending the population of models to include a realistic DRG neuron geometry would allow new research questions to be considered and could extend calibration of the population, for example models could be calibrated against the known range of conduction velocities in C-fibres.The variability in conductance could also be extended to account for differences in ion channel expression at different parts of the axon and cell body [Calvo et al., 2016;Tsantoulas and McMahon, 2014].

Data availability
Zenodo: Data supporting "Changes in human dorsal root ganglion neuron excitability from modulating Nav 1.8 conductance are non-linear and depend on the conductances of the delayed rectifier and M-type potassium currents: a simulation study".https://doi.org/10.5281/zenodo.5512414 This project contains the following underlying data: • drg-pom-master.zip(Source code, models and examples for this study and drgpom software.) • Figure_1_data.csv (Voltage trace data underlying Figure 1A.)

Yes
Is the rationale for developing the new method (or application) clearly explained?Yes Is the description of the method technically sound?Yes Are sufficient details provided to allow replication of the method development and its use by others?Yes If any results are presented, are all the source data underlying the results available to ensure full reproducibility?Yes Are the conclusions about the method and its performance adequately supported by the findings presented in the article?Yes Competing Interests: No competing interests were disclosed.
Reviewer Expertise: microelectrode arrays, electrophysiology, nociceptors I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Markos Xenakis
Institute of Neurophysiology, RWTH Uniklinik, RWTH Aachen University, Aachen, NRW, Germany In the present work, the authors, Oliver J. Britton, and Blanca Rodriguez, make a decisive step toward whole-hDRG-neuron modeling and simulation.They are motivated by the need for more realistic insilico models for studying the transmission of nociceptive stimuli.They systematically investigate how inter-neuronal ion channel density heterogeneity contributes to the shape of hDRG action potential.This is done in a computationally efficient manner without sacrificing biophysics realism.Namely, the authors introduce a pool of hDRG models where model-to-model variability of ion channel conductances but not ion channel kinetics is implemented.Their modeling procedures are documented pedagogically.They emphasize the role of the NaV1.8 channel in shaping action potential properties.Specifically, they report on a balancing mechanism regulating neuron firing rate concerning an inflection point of the NaV1.8 conductance.Driving NaV1.8 conductance above its inflection point value decreases the neuron firing rate.The location of the inflection point is reported to depend on the delayed rectifier potassium conductance.The authors conclude that synergies among the NaV1.8 and other channels, especially potassium, play an essential role in modulating the NaV1.8 channel in nociceptive neurons.
The study paves the way for the replacement, refinement, and reduction of animals in research.
The presented model shows good potential in reproducing hDRG electrophysiological signatures.
It is constructed comprehensively and transparently based on previous works and experimental data.What I missed, however, is a knowledge gain in neuropathic pain disease.Remarks and comments for improvement are attached below.

Major comments:
Methodological aspects: The authors do a great job documenting their procedures and sharing step-by-step guidelines on installing and running the software.However, all this occupies too much space in the manuscript, which could be used for other, more relevant purposes, such as explaining more in-depth how presented findings can help better understand neurophysiological aspects of nociception.I hence suggest that text and tables from p. 7 "Installing drgpom and compiling models .. " to p. 13 ".. .This example can be used as a template to design your own simulation studies." to be placed elsewhere.For example, the authors could consider the option of moving this material to the corresponding Zenodo and/or GitHub repository.

Results:
On a few occasions, it is unclear whether findings reflect specific parameter choices or actual neurophysiological traits.A better explanation of the relationship between the former (parameter space) and the latter (neurophysiology) could help.Specifically: "For Gnav 1.8, this is because models with low GNav 1.8 have peak voltages that are below the experimental range, and models with high GNav 1.8 have threshold voltages that are below the experimental range.For GKA, no model with high GKA had a combination of threshold voltage, peak voltage, and after-hyperpolarisation time constant that was simultaneously within the experimental ranges of these three AP biomarkers."-What is the neurophysiological rationale in support of these simulation results?
"This correlation is caused by calibrating on a combination of three biomarkers: threshold voltage, peak voltage and AP full width parameters, as removing all of these AP biomarkers from the calibration process also removes the correlation between GNav 1.8 and GKdr and removing any one of them reduces the correlation coefficient.This indicates that in our model the balance of these two currents is important throughout the action potential, from upstroke (threshold, peak voltage) to repolarisation (AP full width)."-What is the principle underlying reported correlation?Is it a computational artifact, or does it reflect a biophysically interesting phenomenon?
"The relative proportions of these currents are however highly variable from model to model.Similarly, the main currents during the upstroke are INav 1.7 (red) and INav 1.8 (blue), but some models have a large contribution by Nav 1.7 and in some models the contribution is very small, while the contribution of INav 1.8 is uniformly large."-Could the authors clarify how NaV1.7 and NaV1.8 contribute differently to the current generation?
Conclusions and Discussion: Based on the abstract and keywords, the reader expects to learn something about neuropathic pain disease.After reading the manuscript, however, this expectation still needs to be met.The role of the NaV1.7, which has emerged as a key regulator of neuropathic pain disease, could also be addressed, especially concerning the pivotal role NaV1.8 is suggested to have here.
Could the authors elaborate on how presented findings and modeling procedures enhance our understanding of neuropathic pain?

Minor comments:
On the first page, in the Methods description section, a full stop is missing: ".. in ionic conductances and action potential morphology" should read ".. in ionic conductances and action potential morphology Are the conclusions about the method and its performance adequately supported by the findings presented in the article?Partly Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Computational Neurophysiology I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
This study uses computational modelling to investigate how the human dorsal root ganglion (DRG) neuron responds when the Nav 1.8 sodium channel is selectively modulated.The study is motivated by the role of the DRG neuron in nociception (pain) and the role of the Nav 1.8 channel in modulating the excitability of that neuron.A population of models is used to investigate phenotypically diverse DRG neurons drawn from the human population.They find that the excitability of the neuron depends not only on Nav 1.8 but also on potassium ion channels which appear to be naturally co-regulated with Nav 1.8.
The proposed model offers practical benefits in reducing the use of animals in pain research, which I applaud.The population of models approach is useful for addressing the issue of natural variation too.

Methods
A substantial part of the manuscript is dedicated to instructions for installing and running the software.Documenting the software is important, but those issues are separate from the scientific question that is being investigated (about ion channels).So I urge the authors to move all of the software-related text (pages 7-13) into a separate software manual or tutorial paper.Software documentation is best shipped with the source code.

Discussion
Upon reading the discussion, I felt no wiser about how co-regulated INav 1.8 and potassium currents might inform better treatments for patients with pathological pain.The discussion of the M-type potassium current seemed to hint that it might it be possible to treat a defective Nav 1.8 ion channel by targeting the potassium channels.I would have like to see that sort of thing explored in greater depth.
Minor Comments pg 3. Suggested rewording: "For example, a gain-of-function mutation in Nav 1.It would be nice to plot the stimulus protocols to spare the reader from looking it up.It is also worth describing how the randomized models were equilibrated prior to being analyzed.Some AP models can take hundreds of beats to equilibrate to a new parameter regime.
p14. "This correlation is caused by calibrating on …" I'm not sure this is the best wording.Instead of saying the correlation is caused by the calibration (which is to say that it is caused by your method), I suggest saying that the correlation reveals the normal relationship between ion channels (as caused by nature).
Figure 3A.The a-h labels on the abscissa are slightly confusing.Better to label them as GNav17, GNav18, etc.
Figure 4A is unneccesary.You could omit it altogether.The idea is clear enough from panels B and C.That said, the effect of Nav 1.8 scaling could be made even clearer by plotting it for only one stimulus amplitude, say 4000pA.
p18. "We hypothesised that one reason for the decrease in firing rate …".
Illustrating the idea with some example AP trains would be helpful.

p25. "To estimate the number of animals currently …"
There is no need to explain the details of the pubmed search, just summarize it.For example: "To estimate the number of animals that could be potentially be replaced with our proposed DRG model, we searched pubmed for all papers in the past five years that described pain research using DRG neurons in rodents.We estimated that 21,000 animals were used annually for this purpose.These estimates do not include animal use in industry".

Are a suitable application and appropriate end-users identified? Partly
Are the 3Rs implications of the work described accurately?Partly Is the rationale for developing the new method (or application) clearly explained?Yes

Is the description of the method technically sound? Partly
Are sufficient details provided to allow replication of the method development and its use by others?

Yes
If any results are presented, are all the source data underlying the results available to ensure full reproducibility?Yes Are the conclusions about the method and its performance adequately supported by the findings presented in the article?Partly Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Computational cardiology I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above. Reviewer

Summary:
This manuscript by Britton O.J and Rodriguez B. represents a comprehensive method article introducing an open-source software designed for simulating human dorsal root ganglion (hDRG) neuron electrophysiology.Initially, a baseline model with eight transmembrane currents representing sodium and potassium channels, was established, followed by the generation of a pool of candidate models using a population approach.This aimed to align them with specific experimental AP biomarker ranges.In addition, the research focused on the Na V 1.8 channel conductance to elucidate its influence on neuronal excitability in relation to other ion channels in nociceptive neurons.Notably, modifications in the Na V 1.8 conductance revealed an interesting behavior: a moderate increase heightened the firing rate, but an excessive increase reduced it.
The study also underscored the importance of delayed rectifiers and M-type potassium conductances in determining neuronal excitability.In particular, perturbations in the delayed rectifier potassium conductance affected the Na V 1.8's impact on the firing rate.The authors performed a detailed simulation analysis on the variation of those as well as other conductances on the action potential firing rates revealing interesting trends, which can be tested experimentally and, in the future, also used for the development of new specific ion channel modulators.

Strengths:
The manuscript is well written.The hypothesis is well-defined and the rationale for the study is clear, giving readers a good idea of what to expect in the subsequent sections.In addition, the methodology employed for capturing inter-neuronal variability is noteworthy.
The model accounts for variations in ionic conductances and action potential morphology, making it more comprehensive and likely more representative of biological reality.
A major strength of the paper is the sharing of their software openly.This move makes the research more trustworthy and lets others build on it.Open tools can also help new researchers learn.Instructions for the use of software are very clear.

Areas of potential improvement:
By focusing solely on the soma, the model may not capture the complete electrophysiological behavior of the entire neuron, especially any interactions between the soma and an axon.This is perhaps beyond the scope of the study and can be acknowledged as a study limitation.
The simplifying assumption that ion channel structures do not vary between different DRG neurons might not hold true for all neurons and conditions.There is ion channel structural variability, which may influence their open probabilities, gating time constants and other functional parameters.This assumption could impact the overall accuracy and applicability of the model and can be also mentioned as a limitation.
I would also consider moving sections of the Methods, which discuss software installation, running etc. towards the end of the paper as an Appendix, in order not to break the flow of the paper.
Upon a detailed evaluation of your manuscript, I appreciate the depth of the research presented.I believe that with certain refinements in the figures, Results, and Discussion sections, the manuscript will greatly benefit in its readability and overall impact.
Here are my specific recommendations: Figure 1.The current name on the right got cut.Analysis: It remains ambiguous if these models were run until they achieved a particular steady-state and then fitted it to the mean values of each parameter based on the data recorded from ex vivo hDRG neurons. 1.
For Figure 2, utilizing smaller data points and a more distinguishable color palette could enhance the clarity of the distribution of conductances across the population.

2.
Figure 3, The "currentscape" visualisation technique implemented in Figure 3 offers a 3. unique insight into the conductances and relative current magnitudes across six models.To further its effectiveness: Opt for color schemes that offer greater contrast to discern the various currents more distinctly.Incorporate a dynamic repetition of action potentials to enrich presentation.Figure 9 Colors for GKM/GKdr 0 and 2 as well as for 0.5 and 1.5 are hard to distinguish.8.
While the procedure involving the stimulation of different levels of block and enhancement of Na V 1.8 conductance is mentioned, a clearer description of the methodology could be helpful.This includes specifics about how 'different levels' were defined and chosen.Currently it is not clear whether they (as well as other conductances) were selected solely based on physiological interneuron variability, action of a pharmaceutical or another modulator or both.

9.
The Na V 1.8''s influence on hDRG neuron excitability analysis could benefit from a section discussing potential future investigations or implications of these findings on possible therapeutic interventions, given the evident role of Na V 1.8 in neuronal firing.

10.
The comparison and usage of previous modeling studies of DRG neurons, particularly from Tigerholm et al. and Choi and Waxman, adds depth to the study.It would be advantageous to highlight any adaptations or modifications made to these referenced studies' methodologies and the reasoning behind these changes.

11.
Throughout the paper: It would be best not to use space in the "GNav 1.8" and similar variables.Ideally one would also use subscript for channel name, i.e., "G_{Nav1.8}".Ideally, voltage-gated ion channel names should include subscript capital "V" and not "v" as well. 12.
Please briefly explain the word "rheobase", which is specific to the ion channel electrophysiology field and may not be known to a general reader.
There is no link for "Choi and Waxman, 2011" paper.15.
"simulating changes in multiple conductances at multiple different levels" -Do you mean levels associated with certain changes of conductances from the base model or something else? 18.
"For GNav 1.8, this is because models with low GNav 1.8 have peak voltages that are below the experimental range, and models with high GNav 1.8 have threshold voltages that are below the experimental range."-please check if this is correct as stated, as it seems that both low and high "GNav 1.8 " have the same effect.

19.
"combination of three biomarkers: threshold voltage, peak voltage and AP full width parameters".Are these parameters?I thought they were properties or simulation outputs.

20.
""That GKdr is actually larger in the rapidly firing group suggests it is not just a lack of total 21.potassium current conductance, and a corresponding lack of outward current to oppose repetitive firing, that causes these models to fire so rapidly."This sentence is confusing, please clarify if possible.Maybe, the comma after "conductance" causes this confusion."Mean animals used for DRG studies in the 16/40 papers where number of animals used was reported or could be inferred from figures was 103".Please reword this sentence as it is confusing from the start.What does "mean animals" mean?22.

Conclusions:
The authors have done an excellent job in pointing out areas for potential improvements in their model.Adding more details about calcium-related currents and pumps and transporters will make the model even better.The authors also did not mention chloride channels, which are known to be present in DRG neurons as well (see e.g., Wilke et al Front Neurosci. 2020;14: 287).Also, showing the full structure of hDRG neurons, such as the connections between axons and the main cell body (soma), will make the model more realistic.I look forward to seeing how the authors will refine their work in future updates and how they will extend their population models to capture hDRG neuron variability in human datasets available in the future.

Are a suitable application and appropriate end-users identified? Yes
Are the 3Rs implications of the work described accurately?Yes Is the rationale for developing the new method (or application) clearly explained?Yes Is the description of the method technically sound?Yes Reviewer Expertise: Ion channel and G protein-coupled receptor molecular modeling and simulations.
We confirm that we have read this submission and believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however we have significant reservations, as outlined above.
The benefits of publishing with F1000Research: Your article is published within days, with no editorial bias • You can publish traditional articles, null/negative results, case reports, data notes and more • The peer review process is transparent and collaborative • Your article is indexed in PubMed after passing peer review • Dedicated customer support at every stage • For pre-submission enquiries, contact research@f1000.com

"
Time taken on 4 cores = 1234 seconds.Current population saved to example_population.pkl.This simulation has finished running." 4.You can now open your population results by following the examples in the example Jupyter notebook.Toopen the example notebook follow the instructions in the "Using drgpom to run and analyse simulations" section above.
model_details = { 'mechanisms':{}} Each ionic current model to be included in the model, and each parameter to be varied within that ion channel model is defined as follows: model_details[ 'mechanisms'][ 'nrn_mech_name'] = { 'param_name':'param_name_in_nrn'} 'nrn_mech_name' is the name of the ion channel model in NEURON.To find this name open the.mod file for the relevant current model in the Models directory of the repository and look for the line beginning SUFFIX.This gives the name of the mechanism in NEURON.

Figure 1 .
Figure 1.Experimentally-calibrated population of 328 human DRG neuron AP models.A) Overlaid voltage traces from every model in the population stimulated with a step stimulus at each model's rheobase.Each model was calibrated against hDRG neuron AP biomarker ranges calculated from data in Davidson et al. [2014].B) Model schematic.The model approximates the hDRG neuron soma as a single compartment Hodgkin-Huxley model with 8 ionic currents.Arrows indicate the usual direction of current flow for each current.Each model in the population has the same 8 ionic currents but different channel densities of each to simulate inter-neuronal variability in ionic conductances.

Figure 2 .
Figure 2. Scatter plots and histograms of the distribution of conductances in the population.Histograms on the diagonal panels show the distributions of individual conductances.Scatter plots show distributions for all pairs of conductances.Conductances are given in terms of a dimensionless scaling factor applied to the baseline conductance values (which are provided in the supplementary material).

Figure 3 .
Figure 3. Conductances and ionic currentscapes for six different models from the population of models.A: Conductance scaling factors for six models selected to show a wide variety of different conductances and ionic current profiles.B-G: Currentscapes [Alonso and Marder, 2019]  showing the fractions of each ionic current active during a simulated step stimulus.Each plot shows, from top to bottom, the voltage trace, the total outward current, the fraction each outward current contributes to the total outward current during the simulation, the fraction each inward current contributes to the total inward current during the simulation, and the total inward current.

Figure 4 .
Figure 4. Effects of co-varying Nav 1.8 conductance and stimulus amplitude on firing rate and AP half width.A) Simulation study design showing the grid search design used over stimulus amplitude and GNav1.8.Each dot represents a simulation on all 328 models in the population with the given Nav 1.8 conductance scaling factor applied to all models, and an 800 ms step stimulus current applied with the given amplitude.B) Mean firing rate of the population for each of the different simulation conditions.Each square represents a simulation of all models in the population with a different pair of values for GNav 1.8 scaling factor and stimulus amplitude.C) Mean AP half width of the population for each of the different simulation conditions.

Figure 5 .
Figure 5. Relationship between AP half width and firing rate for individual models in the population of models under different GNav 1.8 scaling factors.Blue dots show results from individual models in different simulations from the simulation design shown in Figure 4A.Red lines show mean values of plotted results.Descending panels show results from simulations with progressively larger GNav 1.8 scaling factors.

Figure 6 .
Figure 6.GKdr and GKM are key determinants for whether a model fires rapidly under many different combinations of GNav 1.8 and stimulus amplitude.A) Conductance distributions of two subpopulations in the population of models.Blue histograms show conductance distribution for the top 25% of models ranked by the number of simulations in Figure 4 in which they fired at a very rapid rate (> 20 AP/s).Red histograms show the other 75% of models.Histograms of the two subpopulations are normalised as the number of models in each subpopulation were not equal.B) Results of varying GNav 1.8, stimulus amplitude and one of GKdr (top row) or GKM (bottom row).Each square in a plot shows the mean firing rate across the population of models for one simulation.

Figure 7 .
Figure 7. Average firing rate of the population of models when non-Nav 1.8 conductances are varied.Each row of plots shows the results of changing one non-Nav 1.8 conductance to 0%, 50%, 100%, 150% and 200% of each model in the population's original value and rerunning the simulation study design shown in Figure 4A.Only variation in GKdr and GKM show substantial effects on firing rate.Within each plot each square represents the mean number of AP/s fired during a step stimulus averaged across all models in the population.

Figure 9 .
Figure 9. Mean firing rate of the population of models averaged over all scaling factors of GNav 1.8 simulated in Figure 8. Left) Each line shows a different GKM scaling factor, with GKdr scaling factor varying on the x-axis.Right) Each line shows a different GKdr scaling factor, with GKM scaling factor varying on the x-axis.

•
Figure_2_data.csv (Conductance scaling factors for each model in the population of models, underlying Figure 2.) • Figure_3A_data.csv(Conductance scaling factors for the models in Figure 3.) • Figure_3B_model_[1/2/3/4/5/6]_data.csv (Voltage and current traces underlying Figure 3B.) • Figure_4A_data.csv(Simulation parameter data underlying Figure 4A.) • Figure_4BC_data.csv(Mean number of APs; AP half width data; and simulation parameters underlying Figure 4B and 4C.) • Figure_5_data.csv (Half width and AP firing rate data for individual models underlying Figure 5.) • Figure_6A_data.csv(Conductance and model category data underlying Figure 6A.) • Figure_6B_GKdr.csv(Mean number of APs and simulation parameters underlying the top row of Figure 6B.) • Figure_6B_GKM.csv (Mean number of APs and simulation parameters underlying the bottom row of Figure 6B.) • Figure_7_data.csv (Amplitudes, conductance scaling factors and firing rates underlying Figure 7.) • Figure_8_data.csv (Mean number of APs and simulation parameters underlying Figure 8.) • Figure_9_left_panel_data.csv (Mean firing rate and GKdr and GKM scaling factors underlying left panel of Figure 9.) • Figure_9_right_panel_data.csv (Mean firing rate and GKdr and GKM scaling factors underlying right panel of Figure 9.) • Equations_and_conductance_densities.docx (Equations and conductance densities for the baseline model.)Data are available under the terms of the Creative Commons Attribution 4.0 International license.
Report 13 December 2023 https://doi.org/10.5256/f1000research.78322.r208690© 2023 Vorobyov I et al.This is an open access peer review report distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.Dr.Gonzalo Hernandez HernandezDepartment of Physiology and Membrane, University of California Davis, Davis, California, USAIgor VorobyovUniversity of California Davis, Davis, CA, USA I have reviewed the manuscript in detail and recommend it for approval with specific reservations.While the paper demonstrates considerable academic merit, there are minor modifications that, if addressed, could enhance the overall strength and impact of the article.

Figure 4 .
I am not sure if panel A is needed, as the same information is already present in panels B and C, as the same range and individual data points as squares are shown there.Although, I do see that Fig.4Ais referenced multiple times for other data as well.4.

Figure 5 .
Figure 5. top panel.No action potential (AP) firing is visible.Maybe, an inset with zoomed in data on this panel would be helpful.5.

Figure 6 .
Figure 6.What is the difference between dark-and light-red plots in panel A? 6. Figures 7 & 8: Both figures summarise the principal trends recognized in the study.To augment the easiness of their interpretation: Enhance the legibility of the x and y axes, the axes are hard to read as the figure contains multiple panels.7.
Are sufficient details provided to allow replication of the method development and its use by others?YesIf any results are presented, are all the source data underlying the results available to ensure full reproducibility?Yes Are the conclusions about the method and its performance adequately supported by the findings presented in the article?Yes Competing Interests: No competing interests were disclosed.
In either case, the model files that need to be compiled are located in the models directory in your installed drgpom folder.If you have installed Anaconda, these will be located in your Anaconda install directory within the Lib/site-packages/ drgpom/models subfolder.Windows you can open the Anaconda folder in your start menu and click the Jupyter notebook icon in there.This will open the notebook in your default web browser, from there, navigate to the notebook of examples (examples_notebook.ipynb) in the examples folder of your drgpom installation and click to open it.A Jupyter notebook is composed of multiple cells containing either text or code.Each cells containing code can be run by selecting it and clicking the run button at the top of the screen or pressing Ctrl + Enter.All code and text in a notebook can be edited so we recommend making a copy of the notebook to experiment with. https://www.neuron.yale.edu/neuron/static/docs/nmodl/unix.html

Table 1 .
Parameters for a simulation using an existing population of models

Table 2 .
Parameters for creating and calibrating a population of models.

Table 2 .
Continued To create and calibrate a new population of models use example_population_creation.py as a template and alter the parameters in Table 1 above and Table 2 below as required.Simulation parameters are the same as those given in Table 1 for using an existing population of models.

suitable application and appropriate end-users identified? Yes Are the 3Rs implications of the work described accurately? Yes Is the rationale for developing the new method (or application) clearly explained? Yes Is the description of the method technically sound? Yes Are sufficient details provided to allow replication of the method development and its use by others? Yes If any results are presented, are all the source data underlying the results available to ensure full reproducibility? Yes
As we show later, the balance of these currents is not as critical as the balance of GNav 1.8 and GKdr at rheobase, where only a single AP is fired."Please clarify what "critical" means in this context.Similar approaches to this have been used .. " should read "Similar approaches have been used .." p.24: Tigerholm et al. (2015) and Choi and Waxman (2011) hyperlinks do not work.
○ p. 16: "○ p. 22: ".. almost complete .. " should read ".. almost completely .." ○ p. 24: "○ Are a Personally, I like to see the equations included in the main text, if only in abbreviated form, so that the reader can appreciate what is being modeled.However, that is not always practical.p6."Each model was simulated …" p5. "The equations for the baseline model can be found in the supplementary material".