In this section we introduce a probabilistic model for the signalling pathway from receptor activation over the G-protein cycle to the calcium channel inhibition. We explain the in vitro experiments which were performed to validate the modelling results and estimate the parameter values. Moreover, we motivate our choice of parameter values and specify the numerical approach used for solving the system.

### The reaction network

The biochemical reaction network under consideration consists of the following reactions (see Table 1 for an overview and Fig. 1 for an illustration). A ligand L attaches to a receptor R in the membrane, resulting in a receptor–ligand complex RL (reaction ({mathcal {R}}_1)). This binding process is dependent on the concentration of protons (i.e. pH) in the microenvironment of the receptor. The protonation of both the ligand and certain residues in the receptor are important determinants for receptor activation, likely due to the pH-dependent formation of hydrogen bonds and/or salt bridges between ligand and receptor18,19.

This receptor–ligand complex RL activates a trimeric G-protein complex which leads to exchange of GDP by GTP and subsequent dissociation into (alpha)– and (beta gamma)-subunits (reaction ({mathcal {R}}_2)). These subunits activate different signalling pathways. Along with the hydrolysis of GTP, another reaction partner M (e.g. arrestin) emerges (reaction ({mathcal {R}}_3)), which initiates internalisation of the receptor–ligand complex (reaction ({mathcal {R}}_4)). The (beta gamma)-subunit inhibits a membrane calcium channel by binding to it (reaction ({mathcal {R}}_5)) (In other papers20 this is referred to as switching the calcium channel from the “willing” state to the “reluctant” state.). After dissociation of the (beta gamma)-subunit from the calcium channel, a trimeric G-protein complex is reformed, and the calcium channel is opened (reaction ({mathcal {R}}_6)). The internalised receptor (RL_w) is either recycled to the cell membrane (reaction ({mathcal {R}}_7)) or degraded (reaction ({mathcal {R}}_8)). The reaction partner M can itself be degraded (reaction ({mathcal {R}}_9)). The ligand L can vanish before it binds to the receptor, e.g. by degradation or unspecific binding to other extracellular components (reaction ({mathcal {R}}_{10})), or it is degraded intracellularly (reactions ({mathcal {R}}_7) and ({mathcal {R}}_8)). The constitutive G-protein activation is given by reaction ({mathcal {R}}_{11}), where we simply use ({mathcal {R}}_2) without ligand. That is, we model the constitutive binding reaction as a net reaction consisting of first the switching of GTP and GDP, followed by the dissociation of the (beta gamma)-subunit from the GPCR complex.

Given these reactions, the stochastic dynamics of the system are mathematically modelled by a reaction jump process characterised by the chemical master equation21,22. The state of the system is given by a vector

begin{aligned} {varvec{x}}=(x_L,x_R,x_{RL},ldots ) in {mathbb {N}}_0^{11} end{aligned}

counting the number (x_S) of molecules of the different species (Sin {mathcal {S}}), where ({mathcal {S}}) is the set of species under consideration:

begin{aligned} {mathcal {S}}:=left{ L,R,RL,RL_w,alpha _{GDP}beta gamma , alpha _{GDP}, alpha _{GTP}, beta gamma , M, Ca_{on}, Ca_{off}right} . end{aligned}

For each reaction ({mathcal {R}}_j) there is a stoichiometric vector (varvec{nu }_jin {mathbb {Z}}^{11}) defining the net change in the population state ({varvec{x}}) induced by this reaction. That is, each time that reaction ({mathcal {R}}_j) occurs, this leads to a jump in the system’s state of the form

begin{aligned} {varvec{x}} mapsto {varvec{x}} + varvec{nu }_j. end{aligned}

E.g., the stoichiometric vector (varvec{nu }_1) of reaction ({mathcal {R}}_1) is given by (varvec{nu }_1=(-1,-1,1,0,ldots ,0)). The rates at which the reactions occur are given by propensity functions (f_j:{mathbb {N}}_0^{11} rightarrow [0,infty )), which can be found in the right column of Table 1.

The temporal evolution of the system is described by the Markov jump process (({varvec{X}}(t))_{tge 0}), ({varvec{X}}(t)=(X_S(t))_{Sin {mathcal {S}}}), where (X_S(t)) is the number of molecules of species S at time t. We define the probability (p({varvec{x}},t):={mathbb {P}}({varvec{X}}(t)={varvec{x}}|{varvec{X}}(0)={varvec{x}}_0)) to find the system in state ({varvec{x}}) at time t given some initial state ({varvec{x}}_0). Then, the overall dynamics are characterised by the standard chemical master equation given by

begin{aligned} frac{d}{dt} p({varvec{x}},t) = sum _{j=1}^{11} left[ f_j({varvec{x}}-varvec{nu }_j)p({varvec{x}}-varvec{nu }_j,t) – f_j({varvec{x}})p({varvec{x}},t) right] . end{aligned}

(1)

The reaction rate equation characterising the corresponding deterministic reaction system is given by the ordinary differential equation (ODE)

begin{aligned} frac{d}{dt}{varvec{C}}(t) = sum _{j=1}^{11}f_j({varvec{C}}(t)) varvec{nu }_j end{aligned}

(2)

for concentrations ({varvec{C}}(t)=frac{1}{V}{varvec{X}}(t)), where V is the system’s volume. It is well-known that the volume-rescaled Markov jump process (({varvec{X}}(t)/V)_{0le tle T}) governed by the chemical master equation (1) converges to the solution ({varvec{C}}(t)_{0le tle T}) of the ODE system (2) in the limit of large particle numbers, i.e., for (Vrightarrow infty)22.

#### Stochastic vs deterministic approach

The stochastic approach has several advantages over the deterministic one. At first, ODEs are an approximation assuming that the higher moments are trivially given by powers of the first moment. Stochastic modelling is exact in the sense that it takes into account all higher moments. Furthermore, the stochastic approach is closer to reality because it assumes a finite set (discrete number) of molecules, while ODEs consider concentrations and only work as approximations for large particle numbers. So the stochastic model is better suited for modelling a small compartment like an axon terminal with a small number of MORs and G-proteins. For our analysis, we will consider comparatively small numbers of molecules for all species (concretely, 20 MORs and 40 G-proteins, see Table 3), such that the stochastic approach is indispensable. Last but not least, a stochastic model delivers more information than ODEs. E.g., it enabled us to analyse the variances of the trajectories or the probability distribution of certain variables like the number of ligand–receptor binding events, which will be done in “Isolated impact of pH value”.

In many situations, however, the ODE model provides a valid approximation of the rescaled first moment of the stochastic process, ({varvec{C}}(t)approx {mathbb {E}}({varvec{X}}(t))/V), as it is also the case here. This fact will be exploited in “Parameter estimation” where the less complex ODE model instead of the stochastic one will be used for estimating the reaction rates (k_1,ldots ,k_{10}) based on experimental data.

### Laboratory in vitro experiments

In order to validate our model, we performed laboratory experiments measuring G-protein activation and membrane calcium currents in vitro. To determine initial G-protein activation (as reflected by the exchange rate of GDP for GTP), the [(^{35})S]-GTP(gamma)S binding assay was used. Because these experiments require genetic alteration (by transfection) of cells, we performed these measurements in commonly used human embryonic kidney (HEK293) cells. In addition, we extracted data produced by FRET experiments3. These experiments measure ligand-induced G-protein subunit dissociation (which follows G-protein activation). The FRET experiments were used to fit the reaction rates. To mimic the mechanisms underlying in vivo opioid analgesia, we examined calcium currents in sensory neurons harvested from rodents using a patch-clamp protocol (see Supplementary Information for methodological details). The experimental results are shown in Figs. 2 and 5a below, and described in more detail in “Isolated impact of pH value”.

### Parameter estimation

Our model includes eleven previously unknown parameters (k_1,ldots ,k_{11}). The determination of appropriate values for these parameters included two main steps: a rough selection of values based on literature, followed by more precise standard parameter estimation.

In Ray et al.7, it has been shown by means of molecular dynamics (MD) simulations that the ligand binding affinity varies for different ligands and pH values. The protonation of both the ligand and certain residues in the receptor are important determinants for receptor activation, likely due to the pH-dependent formation of hydrogen bonds and/or salt bridges between ligand and receptor18,19. We used the relative changes of the rate constant as determined in Ray et al.7 and chose the following different values of (k_1) for different ligand/pH combinations: (k_1=1.25times 10^{-2},text {s}^{-1}) for fentanyl and pH 6.5, (k_1=2.5times 10^{-2},text {s}^{-1}) for fentanyl/pH 7.4, (k_1=2.5times 10^{-3},text {s}^{-1}) for NFEPP/pH 6.5, and (k_1=5times 10^{-4},text {s}^{-1}) for NFEPP/pH 7.4. As the in vitro experiments used for the parameter estimation were performed in the absence of radicals, the rate (k_{11}), which is responsible for the constitutive G-protein activation ({mathcal {R}}_{11}) , was set to zero, because constitutive activation probably does not play an important role in healthy tissue.

The parameters (k_2,ldots ,k_{10}) of the other intracellular reactions were assumed to depend only mildly (if at all) on the ligand/pH combination in the cellular environment; for each of these parameters a single value has been chosen, independent of ligand and pH. This is a reasonable assumption because we chose an intracellular pH value of 7.4 based on well-known mechanisms of cellular homeostasis: although transient (several minutes) changes of intracellular pH may occur with tissue acidosis, intracellular buffer systems and ion pumps in the plasma membrane will rapidly restore physiological pH to ensure cell viability23. Since most previous studies examined situations of longer-lasting inflammation (up to several days)3,4,5,6, we look at this situation, as well. Using results from Zamponi et al.24, we started by setting the rate constant (k_5) of the central binding reaction between the (beta gamma)-subunit and the calcium channel to (k_5=5times 10^{-2},text {s}^{-1}) and proceeded to arrange the other values relative to it according to what is known in the literature. Comparing previous work24,25 it can be deduced that ({mathcal {R}}_5) happens at a timescale an order of magnitude shorter than ({mathcal {R}}_2) and ({mathcal {R}}_3) (with ({mathcal {R}}_2) being slightly faster than ({mathcal {R}}_3)), so we chose (k_2) and (k_3) five resp. ten times smaller than (k_5). The rate constants of reactions ({mathcal {R}}_4,{mathcal {R}}_6,{mathcal {R}}_9) were assumed to be of the same magnitude as those of ({mathcal {R}}_2) and ({mathcal {R}}_3) (with ({mathcal {R}}_6) being slightly faster). The recycling and degradation of internalised ligand–receptor complexes are much slower, at a level of minutes (Fig. 1 in Williams et al.26) which leads to comparatively small rate constants (k_7) and (k_8) for the reactions ({mathcal {R}}_7) and ({mathcal {R}}_8) of the internalised receptor. The extracellular decay of ligand due to unspecific binding and other incidents (reaction ({mathcal {R}}_{10})) was set to a value at which it showed a first effect on calcium channel inhibition.

After this initial step of selecting rough parameter values based on available information, we fine-tuned the parameters (k_1,ldots ,k_{10}) via standard parameter estimation techniques using the in vitro experimental FRET data (see Fig. 2 below) that consists of four individual time series for G-protein activation for the four cases fentanyl/pH 6.5, fentanyl/pH 7.4, NFEPP/pH 6.5, and NFEPP/pH 7.4. The experimental data was first pre-processed by determining the offset time (time point at which the respective ligand was added), and the linear scaling transformation that maps the number of undissociated G-proteins from the model to the measured FRET signal. Then the residual distance between the solution of the ODE model and the experimental data was minimized by optimally adapting the parameters (k_2,ldots ,k_{10}), starting from the initial values previously chosen (second column of Table 2 termed “Preselected”). Here, the residual is the mean squared distance between ODE solution and data, summed over all four time series. The minimization was done using standard techniques for parameter estimation27,28, within the framework of the software PREDICI29. The resulting parameters values are shown in Table 2, third column termed “Parameter Estimation”. These values of (k_2,ldots ,k_{10}), optimally adapted to all four cases of ligand/pH combination at once, were fixed.

In a final step, individually for each ligand/pH combination, the parameter (k_1) was fine-tuned by minimizing the residual function for each single time-series for fixed (k_2,ldots ,k_{10}) by changing the respective (k_1). The outcome was:

(3)

where colors shall help to identify the respective curves in Figs. 3, 4 and 5a.

Some of the optimal parameter values exhibit mild deviations from the preselected ones, but no stark contrast to the literature was observed; in fact, closer inspection showed that the mean squared deviation between model-based simulation and experiment was significantly reduced by fine-tuning parameters. The resulting best fit is shown in Fig. 3.

It has been tested that, for these parameter values, the residual distance between data and stochastic model is very close to the residual distance for the ODE model.

### Numerical simulations of the reaction network

Simulations of the stochastic reaction network were performed using Python 3. For each combination of rate constants, 500 Monte Carlo simulations were carried out and the arithmetic mean was calculated in order to estimate the percentage of closed calcium channels plotted in Figs. 4 and  8. The initial state numbers of receptors, G-proteins and calcium channels were chosen at a ratio of 1:2:4 (see Table 3). These numbers are only a rough estimate since the exact stoichiometry of binding events in relation to the number of activated second messengers is currently not fully understood at the experimental level30,31. However, these numbers should suffice to get some first impressions of the properties of the reaction network.

As a time horizon for each simulation 1200 s were chosen. The results are presented in “Results”. In order to find the steady state of the dynamics under pure constitutive activation (i.e., ignoring the ligand-induced activation of receptors given by reaction ({mathcal {R}}_1)), a simulation without ligands (or, equivalently, with (k_1=0)) was run. Given a non-zero rate constant (k_{11}) for constitutive G-protein activation, we determined the long-term average number a of closed calcium channels under these conditions. This long-term mean a (rounded to natural numbers) of closed calcium channels was then used to determine the initial state for the dynamics including receptor activation (see Table 3 for the results). To check for normal distribution of the mean, the 500 runs were divided into batches of 50 and the respective means then tested. Anderson–Darling test indicated normal distribution with (Ple 0.05), so the 95%-confidence interval of the t-distribution is shown in the plots.