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

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, 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.

We developed a baseline model, which was a single model with no variability in ionic conductances. The baseline hDRG model contains eight ionic current models: INav 1.7 [Vasylyev et al., 2014]; INav 1.8 [Han et al., 2015]; INav 1.9 [Huang et al., 2014]; IKdr [Sheets et al., 2007]; an A-type potassium current (IKA) [Sheets et al., 2007]; an M-type potassium current (IKM) [Maingret et al., 2008]; a leak potassium current (IKleak); and a hyperpolarization-activated current (Ih) [Kouranova et al., 2008]. All current models used a Hodgkin–Huxley formulation.

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². 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 after-hyperpolarisation 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 ±1.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

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:

pip install drgpom

Alternatively, download the latest version from the github repository at: https://github.com/oliverbritton/drg-pom/, or the version used for this study from the repository linked in the data availability statement. To install, first unpack the download, then navigate to its directory, open a command window and run the command:

pip install.

Now that drgpom is installed, the ionic current model files need to be compiled by NEURON before being used. To do this in Windows, follow these instructions:

https://www.neuron.yale.edu/neuron/static/docs/nmodl/mswin.html

For Linux follow these instructions:

https://www.neuron.yale.edu/neuron/static/docs/nmodl/unix.html

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.

Using drgpom to run and analyse simulations

Included in drgpom's examples folder is a Jupyter notebook, showing various examples of how to run simulations and save, load and view results using drgpom. Jupyter notebook (https://jupyter.org/) is a tool that allows code, data and visualisations to be integrated together in a single interactive environment that runs in a web browser.

If you have installed Anaconda you can start Jupyter notebook either by opening a command window and typing:

jupyter notebook

or alternatively in 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.

Running an example simulation

To run an example simulation, open a command window in your drgpom directory (if using Anaconda this will be in the Lib/site-packages/drgpom/examples folder in your Anaconda directory) and run the script example_simulation.py in Python by entering the command:

python example_simulation.py

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

Creating and calibrating a new population of models

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:

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.

‘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.py file 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.

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.

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.

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).

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).

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.

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 after-hyperpolarization. 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 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 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.

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.

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 ± 0.32) than the other models group (1.41 ± 0.43). GKdr is larger in the rapidly firing group (1.60 ± 0.44) compared to the other models group (1.28 ± 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.

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.

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.

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.

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.

Figure 8. Average firing rate of the population of models when GKM, GKdr, GNav 1.8 and stimulus amplitude are co-varied.

Each panel shows the results of simulations with a specific pair of GKM and GKdr scaling factors applied. Within each panel each square represents the mean number of AP/s fired during a step stimulus averaged across all models in the population with GKM, GKdr and GNav1.8 scaled by the relevant values. GKM increases in each panel from left to right, while GKdr increases from bottom to top.

Four main trends can be seen in Figure 8:

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).

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.

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:

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:

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, 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 activity-dependent slowing in a model of a DRG neuron axon. As human-specific potassium current and hyperpolarisation-activated 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 cell-derived 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:

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].