Modeling and analysis of the molecular basis of pain in sensory neurons

Abstract

Intracellular calcium dynamics are critical to cellular functions like pain transmission. Extracellular ATP plays an important role in modulating intracellular calcium levels by interacting with the P2 family of surface receptors. In this study, we developed a mechanistic mathematical model of ATP-induced P2 mediated calcium signaling in archetype sensory neurons. The model architecture, which described 90 species connected by 162 interactions, was formulated by aggregating disparate molecular modules from literature. Unlike previous models, only mass action kinetics were used to describe the rate of molecular interactions. Thus, the majority of the 252 unknown model parameters were either association, dissociation or catalytic rate constants. Model parameters were estimated from nine independent data sets taken from multiple laboratories. The training data consisted of both dynamic and steady-state measurements. However, because of the complexity of the calcium network, we were unable to estimate unique model parameters. Instead, we estimated a family or ensemble of probable parameter sets using a multi-objective thermal ensemble method. Each member of the ensemble met an error criterion and was located along or near the optimal trade-off surface between the individual training data sets. The model quantitatively reproduced experimental measurements from dorsal root ganglion neurons as a function of extracellular ATP forcing. Hypothesized architecture linking phosphoinositide regulation with P2X receptor activity explained the inhibition of P2X-mediated current flow by activated metabotropic P2Y receptors. Sensitivity analysis using individual and the whole system outputs suggested which molecular subsystems were most important following P2 activation. Taken together, modeling and analysis of ATP-induced P2 mediated calcium signaling generated qualitative insight into the critical interactions controlling ATP induced calcium dynamics. Understanding these critical interactions may prove useful for the design of the next generation of molecular pain management strategies.

Figure 1. Schematic of calcium signaling network used in this study.

Ca can enter the cytosol via P2X channels, inositol trisphosphate receptors (IP3R) and passive Ca leakage. ATP binding to P2X activates the channel and induces a rapid increase in cytosolic Ca in the presence of extracellular calcium. ATP binding to P2Y receptors activates membrane-bound phospholipase C (PLC) which hydrolyzes phosphatidylinositol-4, 5-bisphophate (PIP2) into inositol 1,4,5-trisphosphate (IP3) and diacylglycerol (DAG). Cytosolic calcium and IP3 binding triggers the opening of IP3R channels and the subsequent release of endogenous Ca from the Endoplasmic Reticulum (ER) into the cytosol. Cytosolic Ca is translocated to the extracellular medium by plasma membrane Ca ATPase (PMCA) pumps, Na/Ca exchangers (NCX) and to the ER by Sarcoplasmic/Endoplasmic Reticulum Ca (SERCA) ATPase pumps. Phosphoinositides (PIs) are recycled between the plasma membrane and cytosol by phosphorylation and dephosphorylation events. The specific reactions, kinetic constants and non-zero initial conditions used in this study are given in Table 1 and Table 2, respectively.

Figure 2. Coefficient of Variation (CV) of parameters (reaction rate constants and non-zero initial conditions) in the ensemble.

Thirty-one parameters were constrained with a CV of less than or equal to 0.5 and 108 had a CV of less than one. The minimum CV was 0.18 while the maximum was 6.5.

Figure 3. Comparison of model simulations versus training data.

The dashed lines in each case denote the mean simulated value over the ensemble of model parameters while the shaded regions denote one ensemble standard deviation (N = 123). Experimental data are shown with error bars. In each corner, the fraction of experimental points captured at one and three standard deviations is given. (A,B): Steady state fraction of open IP3R channels as a function of cytosolic Ca (A) and IP3 concentration (B). The experimental data was reproduced from Bezprozvanny et al. and Watras et al. , respectively. (C,D): Time-resolved measurements of PIP (C) and PIP2 (D) levels following GPCR activation in SH-SY5Y cells. The PIP/PIP2 data was reproduced from Willars et al., . (E): ATP-induced transient increase in cytosolic Ca following P2X receptor activation in P2X3-transfected GT1 cells. Experimental data reproduced from He et al., . (F): ATP-induced transient increase in cytosolic Ca following P2Y receptor activation in Neuro2a cells. Experimental data reproduced from Lakshmi et al., . (G): ATP-dose dependent fraction of gated P2X3 channels for control (black) and cells treated with GDP--S (blue) from rat DRG neurons. Experimental data reproduced from Gerevich et al., . (H): UTP-dose dependent increases in peak cytosolic Ca levels in rat DRG neurons. Experimental data was reproduced from Sanada et al., .

Figure 4. Fraction of gated P2X3 channels (top) and PIP2 levels (bottom) as a function of time and P2Y activation following P2X and P2Y activation with 10 M ATP.

The height of each bar denotes the ensemble mean while the error bars denote the standard error computed over the ensemble. A: Directly following the addition of ATP, the fraction of gated P2X3 channels and PIP2 levels are at a maximum despite P2Y activation. (B,C): The levels of gated P2X3 channels and PIP2 at 30 s (B) and 60 s (C) decreased relative to the wild-type. D: Gated P2X3 channels and PIP2 levels 60 s after the cessation of P2Y activation relax to their initial levels.

Figure 5. Predicted time course of total inositol phosphate levels (sum of IPx) versus experimental measurements in SH-SY5Y cells.

The dashed line denotes the mean simulated value over the ensemble of model parameters while the shaded region denotes one ensemble standard deviation (N = 123). Experimental data are shown with error bars. The data was reproduced from Willars et al., where muscarinic receptors (another class of G protein coupled receptor) was activated by carbachol in the human neuroblastoma cell line, SH-SY5Y . Both the simulation and experiment were conducted with saturating levels of agonist. No parameters were adjusted for this comparison.

Figure 6. Sensitivity analysis as a function of model output and activation conditions.

Squares denote rate constants while circles denote initial conditions organized by biological function. The mean values of the sensitivity coefficients calculated over the parameter ensemble are shown. Vertical and horizontal lines denote the top 10% of sensitive parameters or parameter combinations. Parameters in the shaded regions are highly sensitive regardless of conditions. A: Comparison of cytosolic calcium sensitivity for P2X versus P2Y activation (100M ATP). B: Comparison of the sensitivity of gated P2X and IP3R channels for P2X receptor activation (100M ATP). C: Comparison of the sensitivity on PIP2 and Gq.GTP levels when both P2X and P2Y receptors were activated (100M ATP). (D,E,F): Average rank-ordering of parameter sensitivities as a function of receptor activation state.

Figure 7. Multi-objective thermal ensemble algorithm used in this study.