We extended a previously developed and validated model of inflammation-induced sensitization of muscle nociceptors to include intracellular interactions that contribute to hyperalgesic priming. The existing model simulates the effects of inflammatory mediators, such as PGE₂, on neuronal responses via the independent activation of two distinct GPCR-mediated pathways, resulting in the activation of PKC and protein kinase A (PKA) within the neuron (Nagaraja et al., 2023). One of the GPCRs activates adenylyl cyclase, which catalyzes the activation of adenosine 3′,5′-cyclic monophosphate (cAMP). Under normal conditions, cAMP primarily activates PKA. However, after inflammation, cAMP also activates Epac, which in turn can activate PKC, thus providing feedback between the two GPCR-mediated signaling pathways (Hucho and Levine, 2007). To account for this feedback, we extended the model by adding descriptions of Epac kinetics (Figure 1A) and a new ion channel, i.e., TRPV4, based on evidence of its contribution to mechanical sensitivity in nociceptors (Ho et al., 2012; Mickle et al., 2015). The equations used to describe the Epac signaling and TRPV4 kinetics are provided in the Supplementary material.

Our extended model of inflammatory pain consists of 59 ordinary differential equations (ODEs) and 141 parameters. Each equation represents one variable in the model, where a variable represents activation or inactivation factors of 18 transmembrane proteins; the intracellular concentrations of K⁺, Na⁺, Ca²⁺, and inositol trisphosphate (IP₃); the endoplasmic reticulum (ER) Ca²⁺ concentration; the active and inactive subunits of the two GPCRs; three membrane-associated enzymes; the concentrations of 18 intracellular proteins and second messenger molecules; and the membrane potential (V_m). We provide a list of the model variables, their descriptions, and initial values in Supplementary Table 1. Using a lumped Hodgkin–Huxley-type formalism (Hodgkin and Huxley, 1952), we calculated changes in V_m at a given time point based on the changes in the currents of all neuronal transmembrane proteins listed above as follows:

where C_m denotes the membrane capacitance and I represents the current through the different transmembrane proteins (described by the subscripts). The model simulates the time course of changes in the variables over 96 h in response to a mechanical force applied before and after two subsequent exposures to an inflammatory mediator (Figure 1B). We used 141 parameters to describe all the modeled mechanisms, and provide a list of the model parameter numbers (used to keep track of the parameters in our simulations), names, values, descriptions, units, and sources of the computational or experimental study from which we adapted or derived their values in Supplementary Table 2. We modified a subset of the model parameters (designated as “modified” in Supplementary Table 2) to match the inflammation-induced changes in AP firing magnitude from literature data (see “Model calibration and validation” section below). In all simulations, we maintained the extracellular concentrations of Na⁺, K⁺, and Ca²⁺ as well as the volume of the nociceptor nerve ending and its C_m at constant values (Nagaraja et al., 2023). We provide the ODEs and other equations describing all the modeled mechanisms, as well as the Nernst potentials and ionic balances for the intracellular concentrations of Na⁺, K⁺, and Ca²⁺ in the Supplementary material.

To drive the model, we provided as input a rectangular pulse with a mechanical force of either 5 or 100 mN (Figure 1B). We selected these values to represent inputs from both the innocuous and pain-inducing ranges of mechanical forces for mouse nociceptors (Cavanaugh et al., 2009; Deuis et al., 2017). We applied three pulses for a period of 60 s each, starting with the first pulse at the 1-h simulation time point. Next, at both the 47-h and 95-h simulation time points, we provided as input a 30-min rectangular pulse of an inflammatory mediator at a concentration of 100 nM, immediately followed by another 60-s rectangular pulse of 5 or 100 mN mechanical force. We chose the 30-min period to evaluate the inflammation-induced sensitization of the neuron’s response to mechanical forces based on literature data, which showed that the threshold reduction due to inflammation peaked between 30 min and 1 h in rat and mouse neurons (Hendrich et al., 2013). At the end of each simulation, our model generated a 96-h time course for each of the 59 model variables.

In all our computational analyses, we focused on the V_m time course, from which we calculated the total number of APs generated following the application of each pulse of mechanical force before and after the addition of the inflammatory mediator. We defined an AP as a V_m spike that crossed over 0 mV and reached a minimum value of +10 mV. We used the MATLAB function FINDPEAKS to identify the APs and to record their heights, widths, as well as the simulation time points at which they occurred. Next, we calculated two fold-change values (FC1 and FC2) for AP firing by dividing the total number of APs generated in response to the application of the force after the first and second inflammatory events (47 h and 95 h; Figure 1B), respectively, by the corresponding value before inflammation. The fold changes in the magnitude of AP firing following the two inflammatory events calculated using the nominal parameter set represented the baseline inflammation-induced sensitization. In addition to input forces of 5 and 100 mN, we calculated FC1 and FC2 values for force inputs of 10, 20, and 50 mN. We performed all computations in the software suite MATLAB R2015b (MathWorks, Natick, MA, USA) and solved the model equations using the MATLAB solver ODE15s with default tolerance levels.

We performed a sensitivity analysis to identify proteins and intracellular molecules and their associated signaling processes whose alterations were key for the development of hyperalgesic priming in the neuron and regulation of subsequent inflammation-induced changes in neuronal AP firing. First, we assessed the model’s robustness by performing a local sensitivity analysis (LSA), where we varied the model parameters near their nominal values (± 1%) as previously described (Nagaraja et al., 2014). Second, we performed a global sensitivity analysis (GSA) to account for the known heterogeneity in the expression of the various proteins at different nerve endings of muscle nociceptors. This analysis allowed us to account for the known variability in the conductance, activation, and inactivation gating factors of the same membrane proteins under different stimuli (Gold and Gebhart, 2010).

To this end, we simulated 50,000 distinct nociceptors and stimulated them with two inflammatory events followed by a mechanical force. First, we generated 50,000 unique parameter sets by randomly selecting parameter values from a fourfold range (twofold in each direction) around the nominal parameter values. To generate the random parameter sets, we used Latin hypercube sampling (MATLAB function LHSDESIGN) (Nagaraja et al., 2014). Next, we performed two sets of simulations using the 50,000 parameter sets, where we drove each simulation using a mechanical force of either 5 or 100 mN applied once before and twice 30 min after the addition of an inflammatory mediator at the 47- and 95-h simulation time points (Figure 1B). We interrupted and eliminated the time course simulations of V_m that did not reach the 96-h time point within 15 min of computation time (wall-clock) or that required time steps smaller than 1 × 10^–12 s. We used this lack of convergence in the simulations to flag parameter sets that resulted in non-physiological kinetic behavior. Accordingly, we only used the simulations that ran to completion to calculate the fold changes FC1 and FC2 described above in the number of APs fired (in response to either 5 or 100 mN force) after each of the two inflammatory events. We used the calculations of APs, FC1, and FC2 to further eliminate the simulations in each force group where the number of APs fired in response to a mechanical force in the absence or presence of an inflammatory mediator was zero. Finally, we eliminated the simulations in which the AP fold-change values, FC1 and FC2, after each of the two inflammatory events were equal to one. We assumed that these simulations represent “non-nociceptive neurons” because their responses to noxious mechanical forces did not change after two consecutive inflammatory events, i.e., they did not have the ability to sensitize. Using the results from the remaining simulations, we performed two analyses, a partial rank correlation coefficient analysis and a parameter distribution analysis.

For each transmembrane protein identified as key for FC1 and FC2 changes, we performed simulations to mimic specific modifications associated with that protein which are likely to occur following the first inflammatory event that “primes” its response to future inflammatory events. In these simulations, we either reduced or increased the expression of a transmembrane protein at 48 h into the simulation (i.e., after the first inflammatory event). To simulate an increase or decrease in a protein’s expression, we multiplied or divided, respectively, the current in Eq. (1) corresponding to that protein by a factor of 2. First, we simulated the protein modifications described above using a model with the nominal parameter set, which represented an average nociceptive muscle afferent neuron. Then, to verify that we could reproduce the effects of different modifications in a population of neurons, we repeated the simulations for every modification using parameter sets from simulations that were classified as primed and non-primed neurons. As in the GSA, we interrupted the simulations that did not reach the 96-h time point of the V_m time course within 15 min of computation time (wall-clock) to flag parameter sets where a modification resulted in non-physiological kinetic behavior. Of the simulations that converged successfully, we calculated the mean ± 1 standard error (SE) of the AP fold changes (FC1 and FC2) after each inflammatory event for every modification. Finally, we performed a Wilcoxon rank sum test to compare the mean values of each simulation with a modification to the corresponding simulation without any modification.

To ensure that the model accurately captured the effect of multiple inflammatory events on the magnitude of the AP firing response to mechanical forces, we calibrated the model using data obtained from electrophysiological measurements in rat spinal neurons with peripheral terminals in the knee joint (Bar et al., 2004). The calibration procedure resulted in fold changes in the number of APs fired by the neuron in response to an innocuous (i.e., 20 mN; Figure 2A) and a noxious (100 mN; Figure 2B) mechanical force input after two subsequent applications of inflammatory mediators to fall within one SE of the experimental data. We designated the final set of model parameter values obtained after this calibration procedure as the nominal parameter set. Figure 2C shows the simulated APs for the case where we stimulated the neuron by a 20 mN force before (i.e., at the 1-h time point) and at 30 min after the addition of PGE₂ at the 47- and 95-h simulation time points. In addition, based on the results of the LSA (performed using the nominal parameter set), we found that V_m was not overly sensitive (sensitivity indices > 100) to any of the model’s 141 parameters, suggesting that the model was stable and robust to small perturbations (± 1%) of its nominal values.

To validate the model, we simulated the AP response (based on the nominal parameter set) to a mechanical force of 100 mN applied before (i.e., at the 1-h time point) and at 30 min after the administration of 100 nM PGE₂ at the 47- and 95-h simulation time points and determined the number of APs fired after each stimulus (Figure 3A). We first compared the model predictions of the fold changes (FC1 and FC2) with the corresponding experimental data derived from rat gastrocnemius muscle neurons (Figure 3B). In addition, we compared the model predictions of the percentage reduction in the mechanical threshold (i.e., the minimum force required to elicit an AP from a neuron) 30 min after the addition of 100 nM PGE₂ at the two time points in the simulation with the corresponding data in rat gastrocnemius muscle neurons (Hendrich et al., 2013). In the model, before the addition of the inflammatory mediator, the minimum force needed to elicit an AP was 1 mN. In both comparisons, our predictions were within ± 1.96 SE of the validation experimental data for all the comparisons, indicating that we cannot differentiate between the average results of the experiment and the model predictions (Figure 3C, error bars show 1 SE of experimental data).

Upon validation, we used the model to establish the baseline inflammation-induced increase in AP firing after two inflammatory events. The AP fold-change values for the individual forces after the first inflammatory event (FC1) were 11 for 5 mN, 10 for 10 mN, 8 for 20 mN, 2.1 for 50 mN, and 1.8 for 100 mN (Figure 3D, solid bars). The AP fold-change values for the individual forces after the second inflammatory event (FC2) were 66 for 5 mN, 43 for 10 mN, 21 for 20 mN, 1.75 for 50 mN, and 3.5 for 100 mN (Figure 3D, open bars). The fold-change values after either inflammatory event decreased as the intensity of the mechanical force stimulus increased.

For the next set of simulations, we selected 5 and 100 mN to represent a mechanical force input from the innocuous and pain-inducing range, respectively, for mouse nociceptors. In addition to simulating two inflammatory events using the nominal parameter set, we repeated the simulations using 50,000 virtual neurons. Of the 50,000 simulations performed for both the 5 and the 100 mN mechanical force inputs, 48,588 and 48,651, respectively, ran successfully. After two additional filtering procedures (see section “Sensitivity analysis” in “Materials and methods”), we extracted 4,310 and 4,027 simulations from the 5 and 100 mN groups, respectively, to use for further analysis. Lastly, in both groups, we classified the simulations as either “primed” or “non-primed” neurons. For the 5 mN simulation set, we identified 3,191 and 1,119 simulations as primed and non-primed neurons, respectively, and for the 100 mN simulation set, we identified 2,101 and 1,926 simulations as primed and non-primed neurons, respectively. Figures 4A, C shows the number of APs fired before (pink circles) and after the first (green circles) and second (blue circles) inflammatory events for the primed neurons in both the 5 and 100 mN groups, where the insets show the average APs fired (horizontal black bars). For both the 5 and 100 mN force inputs, the average fold changes in the number of APs fired in response to the second inflammatory event compared to those fired in the absence of inflammation (FC2; Figures 4B, D, solid bars) were higher than those for the first inflammatory event (FC1; Figures 4B, D, open bars). Overall, these results show that our simulations were able to generate a population of neurons in which hyperalgesic priming occurs in response to both innocuous and noxious stimuli.

To identify the transmembrane proteins that strongly regulated the AP response (specifically the number of APs generated) following the addition of an inflammatory mediator to a pre-inflamed muscle nociceptor across many different nociceptor-signaling conditions, we used two distinct analyses (PRCC and parameter distribution; see section “Sensitivity analysis” in “Materials and methods”). We used the AP fold-change values (FC1 and FC2) and the respective parameter values for each simulation group to perform the PRCC and parameter distribution analyses. For the simulations in which we used a mechanical force of 5 mN as input, the PRCC analysis results showed that for the primed neurons, the model parameters associated with Kv7.2 yielded high and statistically significant correlations (p < 0.01) with the FC1 values (Figure 5A, solid black bar), and the model parameters associated with proteins NaK and TRPA1 as well as molecular processes of G_αq-coupled receptor phosphorylation and phosphorylation of Nav1.8 and Nav1.7 yielded high and statistically significant correlations (p < 0.01) with the FC2 values (Figure 5B, solid black bars). Moreover, for the non-primed neuron simulations, the PRCC analysis results showed that the model parameters associated with proteins NaK, Piezo2, and Nav1.7 yielded high and statistically significant correlations (p < 0.01) with the FC1 values (Supplementary Figure 1A, solid black bars), and the model parameters associated with proteins TRPA1 and Nav1.7 and G_αq-coupled receptor phosphorylation yielded high and statistically significant correlations (p < 0.01) with the FC2 values (Supplementary Figure 1B, solid black bars).

For the simulations in which we used a mechanical force of 100 mN as input, the PRCC analysis results showed that the model parameters associated with Kv7.2 yielded high and statistically significant correlations (p < 0.01) with the FC1 values (Figure 5C, solid black bars), and the molecular processes of G_αq-coupled receptor phosphorylation and activation, phosphatidylinositol 4,5-bisphosphate (PIP₂) activation, and Nav1.8, and Nav1.7 phosphorylation by PKA and PKC yielded high and statistically significant correlations (p < 0.01) with both the FC1 and FC2 values (Figures 5C, D, solid black bars). For the non-primed neuron simulations, the PRCC analysis results showed that the model parameters associated with TRPA1 yielded high and statistically significant correlations (p < 0.01) with the FC1 values (Supplementary Figure 1C, solid black bar), and the model parameters associated with the activation of proteins Nav1.8 and Nav1.7 yielded high and statistically significant correlations (p < 0.01) with the FC2 values (Supplementary Figure 1D, solid black bars).

In the parameter distribution analyses, we calculated the 141 Bhattacharyya coefficients to determine the overlap between the distributions of the normalized values of the model’s 141 parameters in the primed and non-primed neuron groups for both the 5 and 100 mN simulations. While none of the parameters had a considerably small overlap area between the two groups of 141 parameters, the parameters that demonstrated the five lowest values in the 5 mN group were associated with activation or inactivation of ion channels Nav1.7, Kv1.1, and TRPA1 and the rates of Nav1.8 and Nav1.7 phosphorylation by PKC and PKA post-inflammation. In the 100 mN simulation group, the parameters that demonstrated the five lowest values were associated with activation or inactivation of ion channels TRPA1, TRPV4, and Nav1.7 and the rates of TRPA1 phosphorylation by PKC and PKA post-inflammation. Supplementary Tables 3, 4 provide a list of the 141 Bhattacharyya coefficients for all model parameters for the 5 and 100 mN simulation groups, respectively. Figure 6 shows representative examples of two such parameters from each simulation set. For the parameters representing Nav1.7 activation (in the 5 mN force simulation group) and TRPA1 activation (in the 100 mN force simulation group), a larger percentage of the simulations in the primed neuron group fell in the lower range of their normalized values compared to those of the non-primed group (Figures 6A, C, solid vs. dashed lines), indicating that Nav1.7 and TRPA1 activation might be downregulated in primed neurons. Conversely, for the parameters representing Kv1.1 activation (in the 5 mN force simulation group) and Nav1.7 inactivation (in the 100 mN force simulation group), a larger percentage of the simulations in the primed neuron group fell in the higher range of their normalized values compared to those in the non-primed neuron group (Figures 6B, D, solid vs. dashed lines), indicating that these processes might be upregulated in primed neurons. Finally, we combined the results of both analyses and identified five ion channels (Kv7.2, Kv1.1, Nav1.7, Nav1.8, and TRPA1) and two molecular processes (G_αq-coupled receptor phosphorylation and Nav1.7 and Nav1.8 phosphorylation) whose modifications could potentially regulate hyperalgesic priming of mouse muscle nociceptors.

To quantify the extent to which different inflammation-induced neuroplastic modifications in membrane proteins affect subsequent neuronal response, we performed simulations where we modified each model-identified key protein, one at a time, after the first inflammatory event in the primed and non-primed neuron groups for both the 5 and 100 mN force inputs. Specifically, for the three proteins Nav1.7, Nav1.8, and TRPA1, we simulated the effect of their overexpression, and for Kv1.1 and Kv7.2, we simulated the effect of reducing their expression on AP firing. We used these specific modifications based on our model results and literature evidence (Nicol et al., 1997; Blair and Bean, 2002; Kwan et al., 2009; Hameed, 2019). We also performed simulations where we combined all protein modifications to occur simultaneously. We induced these modifications after the first inflammatory event, i.e., at the 48-h simulation time point. We then compared the mean AP fold changes (FC1 and FC2) due to each of these modifications to the corresponding value in the simulations with no modifications (i.e., the baseline sensitization) (Figure 7).

Of the simulations with protein modifications, we identified 2,438 and 868 that ran successfully in the primed and non-primed neuron groups, respectively, when 5 mN force was used as input. When 100 mN was used as the force input, we identified 1,605 and 1,800 simulations that ran successfully in the primed and non-primed neuron groups, respectively. When either an innocuous (i.e., 5 mN) or a noxious (i.e., 100 mN) force was used as input, in the set of primed neurons (i.e., neurons in which subsequent inflammatory events increased AP firing without any modification), modifying TRPA1 (Figures 7A, B, vertically dashed bar for FC2), Kv1.1 (Figures 7A, B, horizontally dashed bar for FC2), and Kv7.2 (Figures 7A, B, white bar for FC2) expression as well as the combined modification of all four key proteins (Figures 7A, B, gray bar for FC2) significantly increased the average AP fold change after the second inflammatory event compared to the simulations with no modifications (Figures 7A, B, black bar for FC2). Of note, while the simulations captured the extent to which distinct neuroplastic modifications induced by hyperalgesic priming affected the neuron’s response to subsequent inflammatory events, the effect was quantitatively higher (for all modifications) when an innocuous force was used as input compared to a noxious force (Figures 7A, B, FC2 values).

In the set of non-primed neurons (i.e., neurons in which subsequent inflammatory events decreased AP firing), when either an innocuous force (i.e., 5 mN) or a noxious force (i.e., 100 mN) was used as input, all the implemented neuroplastic modifications significantly increased the average AP fold change after the second inflammatory event (Figures 7C, D) compared to the simulations with no modifications (Figures 7C, D, black bar for FC2). With the exception of Kv1.1 and Kv7.2 expression reduction (Figure 7D, open and horizontally striped bars for FC2), as observed in the primed neuron group, the AP fold changes induced by the different protein modifications after exposure to the second inflammatory event were higher for the innocuous force input (Figure 7C) than for the noxious force (Figure 7D). Overall, based on our simulation results, hyperalgesic priming involving Nav1.7 and Kv1.1 neuroplasticity (compared to Nav1.8, Kv7.2, and TRPA1) appeared to have the largest effect on the sensitivity of muscle nociceptors to mechanical forces during multiple inflammatory events.

Because a majority of the pathways activated by distinct inflammatory mediators converge within the neuron to increase the intracellular concentrations of PKC and PKA, various membrane-bound and intracellular molecules and proteins involved in PKC and PKA signaling have been investigated for their roles in hyperalgesic priming (Ferrari et al., 2013b; Gangadharan and Kuner, 2013; Araldi et al., 2016a,b; Fang et al., 2021). In fact, numerous behavioral studies have demonstrated that activation of PKCε can consistently induce hyperalgesic priming in rats after multiple exposures to inflammatory mediators. The studies demonstrated that the generation of hyperalgesic priming in rat nociceptors depended on proteins downstream of PKC activation, such as cytoplasmic polyadenylation element-binding protein (CPEB), a protein responsible for translation of mRNA in the peripheral terminals of nociceptors. Furthermore, inhibiting downstream targets of CPEB (e.g., α calmodulin-dependent protein kinase II) or inhibiting nuclear transcription factors (e.g., via the administration of an oligodeoxynucleotide antisense) attenuated priming in rats (Aley and Levine, 1999; Reichling and Levine, 2009; Ferrari et al., 2013a,b, 2015; Araldi et al., 2016a). While we do not explicitly model molecules downstream of PKC activation that are involved in the mRNA translation of different membrane proteins, our PRCC analysis identified molecules that are precursors to PKC activation, i.e., G_αq-coupled receptors phosphorylation and PIP₂ activation, as key for increased AP firing in sensitized muscle neurons during a second inflammatory event (Figures 5C, D, black bars), further demonstrating the importance of PKC signaling in inflammation-induced neuronal response changes. However, because a behavioral pain response (i.e., paw withdrawal) was used as an indicator of the change in the neuronal response induced by hyperalgesic priming in the above-mentioned experimental studies, we still do not have electrophysiological correlates to this observed reduction in mechanical input thresholds. To develop new therapies that target the peripheral neurons, we need to understand how their response (i.e., AP firing) as well as the expression and function of key membrane proteins change during multiple inflammatory events.

Many studies that do measure electrophysiological responses of nociceptors have focused on the ability of PKC and PKA to evoke a change in the magnitude of neuronal AP firing by modifying the gating properties of key Na⁺ and K⁺ ion channels (already present on the neuronal membrane), such as Nav1.8 and Nav1.7 (Vijayaragavan et al., 2004; Wu et al., 2012; Tanaka et al., 2016), Kv1.1 (Nicol et al., 1997; King et al., 2014; D’Adamo et al., 2020), and TRPA1 (Kwan et al., 2006; Karashima et al., 2009; del Camino et al., 2010), among others, via phosphorylation. Therefore, a number of studies have investigated Nav1.8 and Nav1.7 modulation in inflammatory pain. For example, nerve growth factor (NGF)-induced thermal hyperalgesia was inhibited in Nav1.8 knockout (KO) mice compared to saline-treated mice (Kerr et al., 2001), and Nav1.8 inhibition via a pharmaceutical blocker showed reduced excitability in inflammatory and neuropathic rat models of pain (Payne et al., 2015). In addition, nociceptor-specific Nav1.7 KO mice showed an increased thermal and mechanical pain threshold and a decrease in AP response to various inflammatory mediators, including complete Freund’s adjuvant, PGE₂, and NGF (Nassar et al., 2004). However, all these studies reported changes in neuronal responses due to modification of ion channels after only one inflammatory event. Therefore, we do not know if and for how long the inflammation-induced changes in the expression and function of these channels last, because none of the studies measured the neuronal responses after a second inflammatory event.

Our analysis allowed us to characterize the neuronal AP responses after two inflammatory events in two distinct scenarios, one where we pre-exposed a group of neurons to an innocuous force (5 mN) and another where we pre-exposed them to a noxious force (100 mN), both after the administration of an inflammatory mediator. We found that, while the number of APs fired were, on average, higher for the 100 mN force input compared to 5 mN, the AP fold changes, especially after the second inflammatory event, were higher when the neurons were stimulated by an innocuous (i.e., 5 mN) force input, suggesting that the effect of hyperalgesic priming on the increase in the neurons’ excitability was higher when the initial force inputs were innocuous. Furthermore, our ability to simulate two inflammatory events led to the identification of proteins and processes key for initiating hyperalgesic priming (i.e., after the first inflammation event) and for regulating it. For example, in agreement with previous electrophysiological studies, our analysis identified Nav1.7, TRPA1, and Kv1.1 as key in driving the increase in the magnitude of AP firing after a single inflammatory event (Figures 5A, C). In addition to these three proteins, our analysis also identified Nav1.8 and Kv7.2 as well as molecular processes of phosphorylation of G_αq-coupled receptor, Nav1.8, and Nav1.7 as key for the increase in AP firing after a second inflammatory event (Figures 5B, D), which current experimental studies did not measure. Overall, our computational analyses identified five proteins that were key for the increase in AP firing to mechanical forces after two subsequent inflammatory events.

While previous studies and our computational analysis have highlighted the role of individual ion channels or protein kinases in regulating pain responses during multiple inflammatory events (Nicol et al., 1997; Baker, 2005; Kwan et al., 2009; Gold and Gebhart, 2010; King et al., 2014; Nagaraja et al., 2021), we were able to compare the relative effect of modifying different proteins and kinases, one at a time, on AP firing in the same set of 50,000 simulated neurons, which is not feasible in experimental studies. For example, in a group of 2,438 distinct primed neurons, with no specific modification implemented, the average inflammation-induced increase in AP fold change in response to 5 mN force was 3.1 (SE = 0.18) and 11.1 (SE = 1.03) after the first and second inflammatory events, respectively. In the same neurons, five of the six distinct modifications we implemented (involving the proteins identified as key by our GSA), one at a time, significantly altered the average AP fold change after the second inflammatory event (Figure 7A). Of these five modifications, three increased the AP fold change by greater than 70%. Specifically, a twofold decrease in Kv1.1 expression, a twofold increase of Nav1.7 expression, and a combined modification of Kv7.2, Nav1.7, Nav1.8, and TRPA1 increased the average AP fold change by 104, 70, and 77%, with mean FC2 values of 22.6 (SE = 1.7), 18.5 (SE = 1.7), and 19.4 (SE = 1.2), respectively (Figure 7A, open, diagonally dashed, and gray bars). We saw similar results when we simulated the modifications in another group of 1,605 primed neurons using 100 mN force as input (Figure 7B). We can use this information while designing new strategies to regulate the effect of hyperalgesic priming. For example, since the modifications of both Kv1.1 and Nav1.7 individually had a significant effect on AP fold change compared to other protein modifications, we could test them as potential targets to regulate the degree of hyperalgesic priming and select the one with the least side effects.

We also simulated the neuroplastic modifications in the non-primed neuron group. We assumed that these simulations represented silent nociceptors, or non-responding sensory neurons, and aimed to determine whether there were any specific neuroplastic modifications initiated in these neurons by the first inflammatory event that could activate them to sensitize during a second inflammatory event. For example, a group of 1,800 neurons that initially (with no modifications) showed no increase in mean AP fold change after the first or second inflammatory events [Figure 7D, black bars, FC1 and FC2 of 0.50 (SE = 0.004) and 0.45 (SE = 0.004), respectively] became sensitized (or acquired nociceptive properties) after the second inflammatory event when we decreased Kv7.2 or Kv1.1 expression after the first inflammatory event. [Figure 7D, open and horizontally dashed bars, FC2 of 0.75 (SE = 0.20) and 0.74 (SE = 0.17), respectively]. Thus, in addition to identifying a panel of key proteins and molecules that could be potential targets for regulating the neuronal excitability after inflammation, our analysis also provided a one-to-one comparison of the efficacy of targeting each key protein in the same population of neurons.

Our computational model has several limitations arising from simplifying assumptions required to capture the complex nature of hyperalgesic priming in muscle afferent neurons. First, due to limited availability of electrophysiological data for the effects of multiple inflammatory events on AP firing in mice, we used data from rat neurons to calibrate and validate the model, which might affect model accuracy. Second, because we only represent two specific GPCR signaling pathways activated by a subset of the inflammatory mediators that can be present in an injured tissue, it is possible that proteins and molecules involved in other pathways not currently accounted for might be important for regulating hyperalgesic priming. However, as previously discussed, many inflammatory pathways are known to converge within the neuron, resulting in an increase in the concentrations of PKC and PKA (Gold and Flake, 2005; Gold and Gebhart, 2010; Gangadharan and Kuner, 2013; Pak et al., 2018), whose effects, including Epac-mediated PKC activation, we do currently model. Therefore, we could incorporate the effects of other inflammatory mediators if and when relevant. Third, in our model, we adopted many parameter values from previous computational studies developed to describe neurons from animals other than mice or from physiological tissues other than muscle (Bennett et al., 2005; Lindskog et al., 2006; Mandge and Manchanda, 2018). While we performed a validation procedure to match our computational simulations to experimental data recorded from rat neurons, we did not directly derive parameters from single ion channel current measurements in mouse neurons. This simplification could impact the accuracy of certain model parameters. Fourth, we chose a 48-h period between the two subsequent inflammatory events based on experimental evidence that shows that the reduction in the mechanical threshold due to the addition of an inflammatory mediator typically recovers in ∼24 h (Hendrich et al., 2013). However, behavioral studies have demonstrated that the effect of priming can persist for up to 1 week (Ferrari et al., 2013b,2015; Araldi et al., 2016a,b). Therefore, the mechanisms underlying the priming effect might differ based on the timing of subsequent inflammatory events in the primed neurons, which we do not currently consider but could potentially account for in future analyses. Finally, our hypotheses regarding the contributions of Kv7.2, Kv1.1, Nav1.7, Nav1.8, and TRPA1 and both Nav1.8 and Nav1.7 phosphorylation to the hyperexcitability of muscle nociceptors after two inflammatory events stem solely from simulations. These hypotheses need to be tested using experiments involving electrophysiological (patch clamp or dynamic) recordings in mice, and potentially even using human induced pluripotent stem cells, where we separately modify each protein or process and assess the effect of the modification on AP firing in response to multiple, sequential inflammatory events. Ultimately, there is always the question of translatability of the hyperalgesic mechanisms across species. Since we cannot perform in vivo investigations on human nociceptors, mouse models provide the opportunity to use genetic approaches to investigate the molecular mechanisms of nociceptive signaling. In addition, both rodents and humans are known to have the ability to sensitize following repetitive injuries and similar functional organization of the spinal cord, making them a useful substitute for human nociceptors (Fitzgerald, 2005; Toossi et al., 2021).