ar
X
iv
:0
81
0.
51
49
v1
[
co
nd
-m
at.
sta
t-m
ec
h]
2
8 O
ct
20
08 Can a microscopic stochastic model explain the
emergence of pain cycles in patients?
Francesca Di Patti1,2 and Duccio Fanelli1,2,3
1 CSDC Centro Interdipartimentale per lo Studio delle Dinamiche Complesse, via G.
Sansone 1, 50019 Sesto Fiorentino, Italy
2 INFN, sez. Firenze, via G. Sansone 1, 50019 Sesto Fiorentino, Italy
3 Dip. Energetica, Florence University, via di Santa Marta 3, 50139 Firenze, Italy
E-mail: f.dipatti@gmail.com
E-mail: duccio.fanelli@unifi.it
Abstract. A stochastic model is here introduced to investigate the molecular
mechanisms which trigger the perception of pain. The action of analgesic drug
compounds is discussed in a dynamical context, where the competition with inactive
species is explicitly accounted for. Finite size effects inevitably perturb the mean-
field dynamics: Oscillations in the amount of bound receptors spontaneously manifest,
driven by the noise which is intrinsic to the system under scrutiny. These effects
are investigated both numerically, via stochastic simulations and analytically, through
a large-size expansion. The claim that our findings could provide a consistent
interpretative framework to explain the emergence of cyclic behaviors in response to
analgesic treatments, is substantiated.
Keywords: Special issue, chemical kinetics, stochastic processes, population dynamics
1. Introduction
Pain in animals, including humans, is triggered by the so-called nociceptors, sensory
neurons that react to potentially damaging stimulus. Neurology textbooks [1] reports
on the cascade of successive reactions which are activated by the so called noxious
stimuli: The peripheral terminals of primary sensory neurons launch the signal, which
is then transmitted to the spinal and supraspinal nuclei and eventually yields to the
activation of a matrix of cortical areas that are deputed to the conscious experience of
pain.
More specifically, the stimulus originating from a bodily harming menace can be
directly processed through transduction by the receptors located on the nerve terminals.
Alternatively, an indirect pathway can take over through the activation of transient
receptor potentials on keratinocytes or the release of intermediate molecules (such as
the ATP) which, in turn, act on sensory neurons receptors. In the following we shall
assume the first scenario to hold, and, though certainly important, disregard other
mechanisms that might be simultaneously in play. In other words, we simplistically

Can a microscopic stochastic model explain the emergence of pain cycles in patients? 2
imagine that pain receptors act as effective gates, channeling the route to the involved
cortical circuits.
Analgesic drugs relieve the pain by interfering with the peripheral and central
nervous system. Drug molecules bind in fact their target receptors, and consequently
inhibit the pain perception. To grasp and visualize the essence of the process, one can
hypothesize that the bound chemical element occludes the path, by impeding the signal
transduction through the channel envisioned above.
Analgesic are commonly used in basic research and clinical practice, but their
interaction with nociceptory and normal sensory processing remains to be fully
unraveled. Anesthetics are for instance known to modify the electrical recordings
measured via evoked potentials (EPs) responses [2], a powerful diagnostic tools employed
to monitor and characterize a large variety of central nervous system disorders. EPs are
elicited by a specific stimulus applied to the e.g. pain receptors and consist in recording
the induced electrical brain activity, as detected by localized electrodes placed on the
surface of the head. Furthermore, EPs are also useful in documenting objective response
to pain [3, 4] and can thus prove fundamental to elucidate the molecular processes that
control anesthetic absorption and metabolization.
Different analgesic agents have been shown to produce intriguingly distinct effects
at the level of the EPs [5]. Recorded time series of the solicited electric activity display
in fact remarkably different patterns, which are generically attributed to the chemical
specificity of the anesthetic compound. Qualitatively, large, regular, oscillations of
the electric response manifest, latency and amplitude being peculiar traits, supposedly
related to the molecular characteristic of the administered drug.
Furthermore, cycles in the perception of pain have been also reported which might
be hypothetically driven by similar microscopic processes, the interaction between
the anaesthetic molecules and their targets playing certainly a role of paramount
importance. Clearly, the individual experience of pain is also influenced by psychological
and cultural factors, unfortunately difficult to deconvolve when aiming at resolving the
objective picture.
The issue of developing a unique interpretative framework to account for the
presence of such oscillatory regimes has catalyzed vigorous discussions. The puzzle
of their existence remains however to be fully understood.
Current mathematical models [6] approach the problem via deterministic
paradigms, thus neglecting the crucial role which is certainly played by the noise,
intrinsic to the phenomenon under scrutiny. These aspects become particularly
important when accounting for the presence of diverse chemical species, which populate
the stream flow in a spatially diffusive environment. Different chemical entities may
compete with the drug molecules and occupy the sites located in close vicinity of the
receptors, thus effectively hindering the binding event. Under specific conditions, such
competition sustained by the stochastic component of the dynamics might result in large
temporal oscillations for the amount of bound receptors, a mechanism which could
explain the emergence of macroscopic cycles for the sensation of pain in response to

Can a microscopic stochastic model explain the emergence of pain cycles in patients? 3
medicaments.
In this paper, we shall speculate on the above scenario by putting forward a network
of chemical reactions and performing a system-size expansion through the celebrated
van Kampen theory [7]. This enables us to derive a set of linear equations for the
fluctuations, with coefficients related to the steady-state concentrations predicted from
the first-order theory (i.e. the deterministic rate equations). Solutions are identified
for which the deterministic steady-state occurs via damped oscillations: the inclusion of
second-order fluctuations leads then to the amplification of sustained oscillations. These
conclusions are briefly discussed with reference to the existing medical literature.
2. Description of the model
Within the simplified scenario depicted above, we shall model the chemical interaction
between a large, though finite, number of drug molecules (anesthetic), hereby termed
T , and free receptors RF which represent their binding target. Following a successful
binding event, a molecule of the species T disappears, leaving an empty case, hereafter
labelled E. The population of bound receptor RT is in turn increased by one unit.
These assumptions formally translate into the compact chemical notation:
RF + T
α−→ RT + E (1)
where α stands for the associated reaction rate. The inverse reaction corresponding
to the spontaneous detachment of the bound component may occur ‡ with a certain
probability β§, which motivates the introduction of the dual relation:
RT + E
β−→ RF + T
The anesthetic molecules T surf in a densely packed environment and certainly
experience collisions with several other microscopic chemical entities, H , which populate
the streaming flow. Binary interactions between H and T elements, can occur in
the close vicinity of the receptors (RF ) location, potentially disturbing and eventually
interfering with the binding event. As a result of an hypothetical collision, the active
species T can be ejected by the spatial layer immediately adjacent to the receptor,
leaving behind an empty case E. Beyond this effect, which stems from purely steric
interactions, one has to account for possible chemical transformations, which might
occur when individuals from the H and T species encounter: The active chaser T
‡ We here assume that the free T molecule is still chemically active and can thus potentially chase for
unscreened targets. This working hypothesis can be relaxed leading to conclusions qualitatively similar
to the ones highlighted below.
§ In principle it would be extremely useful to dispose of experimental estimates for the reaction rates,
so to define a realistic range of variability for the free parameters in the model. The most reliable data
concern the so-called (equilibrium) affinity constant for the case of e.g. the Tramadol, an analgesic
agent which belongs to the class of synthetic opioid. Depending on the target receptor (and on the
specificity of the chaser’s molecule) the affinity constant is reported to vary of a large amount which
scans two orders of magnitude (from fraction of unity to hundreds) [8, 9].

Can a microscopic stochastic model explain the emergence of pain cycles in patients? 4
RF + T
α−→ RT + E
RT + E
β−→ RF + T
H + T
γ−→ H + H
H + T σ−→ H + E
RT RF
RF
TH
γ
σH
T
Figure 1. The main reaction schemes are depicted. The squares stand for the
inactive species, while the circles represent the drug molecules. The model is then
complemented with a set of additional reactions (see equations (2)–(5)), which accounts
for the possibility that T and H enter (resp. leave) the region deputed for the
interaction.
can lose its ability to bind the target ‖ and it is thus mapped into an inactive H
molecule. To incorporate these effects into the proposed description we postulate the
following interaction rules, which are loosely inspired to the predator-prey competition
mechanism:
H + T
γ−→ H + H
H + T σ−→ H + E
The values of σ and γ characterize the effectiveness of the interaction, which is in turn
sensitive to the choice of the compound T . The idealized cartoons of figure 1 are aimed
at visualizing the above reaction schemes.
To complete the model we introduce an effective migration, by requiring that the
T and H molecules can enter (resp. leave) the region deputed to the interaction. The
‖ Note that this can happen both due to a mechanical stress or via chemical combination of the colliding
species, see for instance [10] where the plasma protein binding is discussed. For a specific application
relative to the case of the Tramadol refer to [11].

Can a microscopic stochastic model explain the emergence of pain cycles in patients? 5
latter assumption yields to the following set of chemical relations:
T δ1−→ E (2)
E
η1−→ T (3)
H δ2−→ E (4)
E
η2−→ H (5)
The population, namely the ensemble of elements belonging to an homologous species
X, will be labelled in the following with the symbol nX . Notice that the number of
receptors N1 = nRT + nRF and the total amount of molecules (including the empties)
N2 = nT +nH +nE are conserved quantities. This observation enables us to reduce the
complexity of the problem by setting:
nRF = N1 − nRT nE = N2 − nT − nH
In the following we shall use the vectorial notation n = (nT , nRT , nH) to help keeping
the mathematical developments compact.
We are now in a position to define the transition rates T (n
′ |n) from a state n to
a different one n
′
. In our convention initial states are on the right and final states
on the left. As an example, consider equation (1). The probability to pick up a
T constituents follows from simple combinatorics and reads nT/N2, while there is a
probability (N1 −nRT )/N1 of RF being chosen. This results in α(nT/N2)(N1 −nRT )/N1
for this particular transition rate [12]. A complete listing of transition probability
associated to the preceding set of chemical reactions is here enclosed:
T (nT − 1, nRT + 1, nH |n) = α
nT
N2
N1 − nRT
N1
T (nT + 1, nRT − 1, nH |n) = β
N2 − nT − nH
N2
nRT
N1
T (nT − 1, nRT , nH + 1|n) = γ
nH
N2
nT
N2
T (nT − 1, nRT , nH |n) = σ
nH
N2
nT
N2
+ δ1
nT
N2
T (nT + 1, nRT , nH |n) = η1
N2 − nT − nH
N2
T (nT , nRT , nH − 1|n) = δ2
nH
N2
T (nT , nRT , nH + 1|n) = η2
N2 − nT − nH
N2
Having defined the transition rates, the master equation governing the evolution of
the discrete stochastic model takes the form:
d
dt
P (n, t) = T (n|nT + 1, nRT − 1, nH)P (nT + 1, nRT − 1, nH , t)
+ T (n|nT − 1, nRT + 1, nH)P (nT − 1, nRT + 1, nH , t)
+ T (n|nT + 1, nRT , nH − 1)P (nT + 1, nRT , nH − 1, t)
+ T (n|nT + 1, nRT , nH)P (nT + 1, nRT , nH , t)

Can a microscopic stochastic model explain the emergence of pain cycles in patients? 6
+ T (n|nT − 1, nRT , nH)P (nT − 1, nRT , nH , t)
+ T (n|nT , nRT , nH + 1)P (nT , nRT , nH + 1, t)
+ T (n|nT , nRT , nH − 1)P (nT , nRT , nH − 1, t)
−
[
T (nT − 1, nRT + 1, nH |n) + T (nT + 1, nRT − 1, nH |n)
+ T (nT − 1, nRT , nH + 1|n) + T (nT − 1, nRT , nH |n)
+ T (nT + 1, nRT , nH |n) + T (nT , nRT , nH − 1|n)
+ T (nT , nRT , nH + 1|n)
]
P (n, t) (6)
where P (n, t) is the probability to find the system in the state n at time t. In the
next section we shall shortly report about our stochastic simulations, before turning to
develop the analytical framework.
3. Numerical simulations
Based on the chemical equations introduced above, numerical simulations can be
carried on, which respect the intrinsic stochastic nature of the model. To this end we
employ the celebrated Gillespie algorithm [13] which exploits the information encoded
in the reaction scheme to advance the system in time through a sequence of random
increments¶. A randomly selected reaction is forced to occur during each successive
step. It should be emphasized that time increments and associated reactions are chosen
so to recover the exact probability distribution of the stochastic time series. For a more
detail account on the philosophy of the integration recipe, the reader can consult [13].
A typical result is represented in figure 2 where the normalized population of T -
like molecules and bound receptors RT are reported, as function of time. Notice that
large stochastic oscillations are clearly displayed, despite the relatively large number of
simulated molecules. Even more interestingly, the oscillations persist when increasing
the populations amount and only when a very large number of constituents is allowed,
the system eventually sets down to its mean-field solution.
As already remarked in [12, 14, 15], this result is odds with the intuitive believe
that fluctuations can be safely ignored, due to the usual statistical factor 1/
√
N .
Indeed, the observed fluctuations arise from the amplification of the intrinsic noise,
associated to the discrete component of the dynamics and can potentially bear and
extraordinary conceptual relevance in several applications. With reference to the case
at hand, emerging, regular, oscillatory patterns in the RT quota, could potentially
¶ The standard implementation of the Gillespie algorithm is based on a nested sequence of operations
which is here briefly recalled. First, one initializes the system at t = 0, by assigning a number of
molecules to each of the considered species. Then the following iterative scheme runs: (i) Calculate the
transition rates Ti(n′|n) associated to each prescribed reaction i. The quantity T0 =
∑M
i=1 Ti(n
′|n) is
stored; (ii) Extract two random numbers, respectively r1 and r2, from a uniform distribution, which are
used to (a) update the simulation time by a finite amount δt = 1/T0 ln(1/r1) and (b) select the index i
labelling the next reaction which is deputed to occur (i is such that
∑i−1
k=1 Tk < r2T0 6
∑i
k=1 Tk); (iii)
Update the species accordingly and go back to point (i).

Can a microscopic stochastic model explain the emergence of pain cycles in patients? 7
4000 8000 12000 16000
Time
0,14
0,16
0,18
0,2
0,22
0,24
φ T 0 1000 2000 30000
0,1
0,2
4000 8000 12000 16000
Time
0,22
0,24
0,26
0,28
0,3
0,32
0,34
φ R
T
0 1000 2000 3000
0,1
0,2
0,3
0,4
(a) (b)
Figure 2. Drug molecules (panel (a)) and bound receptor (panel (b)) densities as a
function of time. The solid lines refer to the time series from stochastic simulations,
while the dashed lines are calculated from numerical integration of the mean field
equations (7)-(9). The insets represent a zoom of the initial evolution and allow to
appreciate the agreement with the mean–field solution at short times. Parameters
values for both panels are α = 0.008, β = 0.005, γ = 0.3, δ1 = 0.001, δ2 = 0.05,
η1 = 0.001, σ = 0.06, N1 = 5300 and N2 = 20000.
explain the modulations in the pain perception as reported in the literature, higher
values of nRT corresponding, in our interpretative scheme, to a less pronounced pain
level. Alternatively, on a different time scale, the effects here discussed could provide a
consistent microscopic picture to understand the presence of quasi–periodic fluctuations
in the evoked electric activity of laboratory animals under anesthetic treatment.
In the forthcoming sections, we shall gain some analytical insight into the model
results and characterize the specific traits of the observed oscillatory behaviors through
the elegant van Kampen [7] expansion.
4. On the deterministic limit
The deterministic counterpart of the governing master equation is straightforwardly
obtained as follows. Focus on T and observe that, by definition:
〈nT 〉 =
∑
n
nTP (n, t)
Multiplying equation (6) by nT and summing over n returns on the right–hand
side d 〈nT 〉 /dt. Simplifying the left–hand side is somehow more laborious and requires
some effort. Proceeding in a completely analogous fashion for nH and nRT yields to the
following system of coupled differential equation for the ensemble independent variables:
d
dt
〈nT 〉 = − α
〈
nT
N2
N1 − nRT
N1
〉
+ β
〈
N2 − nT − nH
N2
nRT
N1
〉
− (γ + σ)
〈
nH
N2
nT
N2
〉

Can a microscopic stochastic model explain the emergence of pain cycles in patients? 8
− δ1
〈
nT
N2
〉
+ η1
〈
N2 − nT − nH
N2
〉
d
dt
〈nRT 〉 = α
〈
nT
N2
N1 − nRT
N1
〉
− β
〈
N2 − nT − nH
N2
nRT
N1
〉
d
dt
〈nH〉 = γ
〈
nH
N2
nT
N2
〉
+ η2
〈
N2 − nT − nH
N2
〉
− δ2
〈
nH
N2
〉
The mean–field approximation corresponds to ignore the correlations, in the above
rate equations when performing the limit for N1 and N2 large. Introducing φT =
〈nT 〉/N2, φRT = 〈nRT 〉/N1, φH = 〈nH〉/N2, rescaling time as τ = t/N2 and formally
sending N1, N2 → ∞, one finally obtains:
d
dτ
φT = − αφT (1− φRT ) + βφRT (1− φT − φH)− (γ + σ)φHφT (7)
− δ1φT + η1(1− φT − φH)
d
dτ
φRT = c [αφT (1− φRT )− βφRT (1− φT − φH)] (8)
d
dτ
φH = γφHφT + η2(1− φT − φH)− δ2φH (9)
where c = N2/N1.
We shall hereon limit our discussion to the case η2 = 0 which will prove analytically
tractable. The conclusions here demonstrated with reference to the selected case study,
will obviously apply to the more general setting where fresh H molecules are allowed
to enter the interaction region. Investigations on the complete model (η2 6= 0) will be
object of a forthcoming publication.
Two fixed points, respectively labelled φ∗1 and φ
∗
2, are identified:
φ∗1 =
(
η1
δ1 + η1
,
αη1
βδ1 + αη1
, 0
)
φ∗2 =
(
δ2
γ
,
α[δ2(γ + σ) + η1γ]
α[δ2(γ + σ) + η1γ] + β[(γ + σ)(γ − δ2) + δ1γ]
,
η1(γ − δ2)− δ1δ2
δ2(γ + σ) + η1γ
)
In figure 2, the solid lines represent the above equilibrium solution: For such choice
of the parameter, as previously mentioned, the stochastic dynamic displays macroscopic
oscillation for the monitored quantities, the average reference value being correctly
predicted by the mean-field theory. What is the reason for these regular fluctuation
to manifest? Are they resulting from the intrinsic finite size nature of the simulated
medium?
To answer these questions it is crucial to determine the stability matrix associated
with the fixed points, as it will play a central role in investigating the cycling
phenomenon. One can thus re–write the mean–field equations in the compact form
dφk/dτ = fk(φ), where the index k = 1, 2, 3 codes the different species, namely T
(k = 1), RT (k = 2) and H (k = 3). The 3× 3 Jacobian matrix Jij = ∂fi/∂φj |FP (here
FP means “evaluated at the fixed point”) controls the linearized dynamics about the
fixed point and can be cast in the form:
J(φ∗) =
0
B
@
−α(1 − φ∗RT ) − βφ
∗
RT
− (γ + σ)φ∗H − δ1 − η1 αφ
∗
T + β(1 − φ
∗
T − φ
∗
H) −βφ
∗
RT
− (γ + σ)φ∗T − η1
c[α(1 − φ∗RT ) + βφ
∗
RT
] −c[αφ∗T + β(1 − φ
∗
T − φ
∗
H)] cβφ
∗
RT
γφ∗H 0 γφ
∗
T − δ2
1
C
A
(10)

Can a microscopic stochastic model explain the emergence of pain cycles in patients? 9
One can easily show that J(φ∗1) admits two real negative eigenvalues and a third
one, also real, whose sign depends on the choice of the parameters. A stability analysis
for the second equilibrium point proves technically difficult. However, via numerical
inspections, a large region in the parameters’ space is identified, which yields to complex
solutions. In particular, complex eigenvalues of the J(φ∗2) are found having negative
real part. This condition ensures an oscillatory approach to equilibrium, a fundamental
ingredient which is eventually responsible for the large scale modulations observed in the
finite size regime. Tracing the region in space deputed to the aforementioned behavior
falls out of the scope of the present paper and shall be postponed to a forthcoming
contribution together with a detailed characterization of the general case with η2 6= 0.
Following the above, from now on, we shall denote with φ∗ a particular value of φ∗2
for which damped oscillations do occur in the mean–field scenario.
5. Characterizing the fluctuations: The van Kampen expansion
As clearly depicted in figure 2 the innate discreteness of the stochastic medium drives
into the system important effects which cannot be captured within the continuous mean–
field scenario. To shed light into the observed phenomena we can bring into the game
the fluctuations by performing the following explicit replacement:
nT = N2φT (t) +
√
N2ξT nRT = N1φRT (t) +
√
N1ξRT nH = N2φH(t) +
√
N2ξH (11)
where the new continuous variable ξ = (ξT , ξRT ξRH ) replace the discrete quantity n =
(nT , nRT , nH) in the definition the probability distribution, namely P (n, t) = Π(ξ, τ).
Before proceeding, it is worth emphasizing that the 1/
√
N1 (resp. 1/
√
N2) term
holds because of the central-limit theorem: The fluctuations are in fact expected to decay
in a similar fashion and, in the continuous limit, N1, N2 → ∞ the system is entirely
characterized in term of its continuous variables φ = (φT , φRT , φH) as prescribed by
equations (11).
The master equation can be re-written as function of the new variables. The left-
hand side reads:
d
dt
P (n, t) =
1
N2
∂Π
∂τ
− 1√
N2
d
dτ
φT
∂Π
∂ξT
− c
−1
√
N1
d
dτ
φRT
∂Π
∂ξRT
− 1√
N2
d
dτ
φH
∂Π
∂ξH
where use has been made of the fact that the time derivative is taken at constant n.
The right-hand side follows from a straightforward, though lengthy, application of the
large-N van Kampen expansion. The main step of the derivation are reviewed in the
appendix. The interested reader can refer to [7, 16] for a detailed account on the whole
procedure.
At leading order of the expansion we recover the deterministic mean–field equations
(7)–(9), while at next–to–leading order we obtain the linear multivariate Fokker Planck
equation (A.5) that governs the evolution of the fluctuations. The coefficients of
this equation are completely specified as function of the model’s parameters (see the
appendix): In principle, by solving equation (A.5) we are in a position to quantify the

Can a microscopic stochastic model explain the emergence of pain cycles in patients?10
deviation from the ideal mean–field dynamics, via the probability distribution Π. At
present, we aim at understanding the oscillation and, to this end, we invoke a completely
equivalent formulation of the Fokker Planck equation. The problem can be in fact cast
as a set of stochastic differential equations of Langevin type, which take the explicit
form:
dξi
dτ
= Ai(ξ) + ηi(τ) i = 1, .., 3 (12)
where for convenience, as a natural extension of the notation introduce above, we have
now set ξ1 = ξT , ξ2 = ξRT , ξ3 = ξH and Ai is specified in the appendix. The term ηi is
a Gaussian noise with zero mean and with correlation given by
〈ηi(τ)ηj(τ
′
)〉 = Bijδ(τ − τ
′
)
To highlight the existence of a possible oscillatory behavior we perform a Fourier
analysis of equations (12), and obtain:
−iωξ˜i(ω) =
∑
j
Mij ξ˜j(ω) + η˜i(ω)
where the tilde stands for the Fourier transform. Following [15] we can re-write this as:
∑
j
Φij(ω)ξ˜j(ω) = η˜i(ω) (13)
with Φij(ω) = −iωδij −Mij . In addition one gets:
〈η˜i(ω)η˜j(ω
′
)〉 = Bij(2pi)δ(ω + ω
′
)
Solving equation (13) for ξ˜i and computing the power spectrum results in:
Pi(ω) = 〈|ξ˜i(ω)|2〉 =
∑
j
∑
k
Φ−1ij (ω)Bjk(Φ
†
ij)
−1(ω) (14)
where we have used Φ†ij(ω) = Φji(−ω). The power spectrum predicted by equation (14)
is plotted in figure 3, for the same set of parameters as employed in the simulations of
figure 2. A clear peak is detected. Moreover, the theoretical curve interpolates correctly
the numerical profile. These results confirm that the macroscopic oscillations which
manifest in the acquired time-series stem from the noise intrinsic to the system under
investigation. It is our believe that mechanisms similar to the ones here hypothesized,
are potentially in place in the complex human (animal) body environment and could in
principle form the basis of a consistent molecular interpretation for the large collection
of experimental, biomedical observations to which we made reference in the introductory
section.
6. Conclusions
In this paper we propose a discrete dynamical framework which is aimed at shedding
new light onto the complex molecular processes intervening in response to an external
harming stimulus so to trigger the pain sensation. We are in particular interested in

Can a microscopic stochastic model explain the emergence of pain cycles in patients?11
0 0,01 0,02 0,03 0,04 0,05 0,06
ω
0
20
40
60
80
100
P(
ω
)
Figure 3. A plot of the power spectrum P (ω) for the RT time series as a function
of the frequency ω. The noisy line corresponds to the spectrum calculated from 500
runs of the stochastic simulation of the model. The smooth line shows the analytic
prediction from equation (14). For the parameters’ setting refer to the caption of figure
2.
elucidating the crucial interplay between the administered drug molecules, which express
their analgesic function chasing the target receptors, and other chemical elements freely
diffusing in the stream. The latter can substantially reduce the anaesthetic effect,
by hindering the available binding sites. Similarly, drug molecules can be turned
into inactive species following binary encounters. The mechanisms here postulated
are formally coded via chemical reactions and define a consistent stochastic scheme.
Numerical simulations display macroscopic oscillations in the concentration amount:
the number of bound receptors change cyclically in time, a trend which we assume to
induce an analogous modulation for the experienced perception of pain. These findings
are analytically illustrated by developing a large-size expansion which enables us to
predict the existence of a peaked power spectrum. It is important to remark that
the amplification process here discussed stems from the underlying stochasticity, which
excites the resonant frequencies of the system. Oscillations arise hence naturally, driven
by the noise which is intrinsic to the system and without invoking any ad hoc couplings
among the molecular agents participating to the dynamics. Our findings could suggest
the existence of a simple, though general, molecular mechanism responsible for the
emergence of cyclic behaviors in response to analgesic treatments [5]. We shall here
stress that our main conclusions apply also for the more general setting where η2 6= 0.
In particular the peaked power spectrum is also found in this latter case and the region in
the parameters’ space which corresponds to the emergence of the cycles can be partially
identified on the basis of explicit analytic formulae. These findings will be reported in
a forthcoming publication [17].

Can a microscopic stochastic model explain the emergence of pain cycles in patients?12
Appendix
The first technical point of the van Kampen expansion concerns the introduction of the
so–called shift operators, E±1T ,E±1RT ,E
±1
H which obey:
E±1T f(n, t) = f(nT ± 1, nRT , nH)
E±1RT f(n, t) = f(nT , nRT ± 1, nH)
E±1H f(n, t) = f(nT , nRT , nH ± 1) .
The master equation (6) is hence cast in the form:
1
N2
∂Π
∂τ
− 1√
N2
d
dτ
φT
∂Π
∂ξT
− c
−1
√
N1
d
dτ
φRT
∂Π
∂ξRT
− 1√
N2
d
dτ
φH
∂Π
∂ξH
=
+ (E+1T E−1RT − 1)
[
α
(
1− φRT −
1√
N1
ξRT
)(
φT +
1√
N2
ξT
)
Π
]
+ (E−1T E+1RT − 1)
[
β
(
φRT +
1√
N1
ξRT
)(
1− φT −
1√
N2
ξT − φH −
1√
N2
ξH
)
Π
]
+ (E+1T E
−1
H − 1)
[
γ
(
φT +
1√
N2
ξT
)(
φH +
1√
N2
ξH
)
Π
]
+ (E+1T − 1)
[
σ
(
φT +
1√
N2
ξT
)(
φH +
1√
N2
ξH
)
Π + δ1
(
φT +
1√
N2
ξT
)
Π
]
+ (E−1T − 1)
[
η1
(
1− φT −
1√
N2
ξT − φH −
1√
N2
ξH
)
Π
]
+ (E+1H − 1)
[
δ2
(
φH +
1√
N2
ξH
)
Π
]
(A.1)
The advantage of using the shift operators relies in that they admit a simple expansion
in the limit for N1 (resp. N2) large:
E±1T = 1±N
−1/2
2
∂
∂ξT
+
1
2
N−12
∂2
∂ξ2T
± · · · (A.2)
E±1RT = 1±N
−1/2
1
∂
∂ξRT
+
1
2
N−11
∂2
∂ξ2RT
± · · · (A.3)
E±1H = 1±N
−1/2
2
∂
∂ξH
+
1
2
N−12
∂2
∂ξ2H
± · · · (A.4)
Plugging (A.2)-(A.4) into (A.1), after some algebraic manipulation, one recovers at the
leading order the mean–field equations, formally identical to the ones reported above,
see equations (7)-(9). The next–to–leading order result in a Fokker Planck equation for
the probability distribution Π(ξ, t):
∂Π
∂τ
= −
∑
i
∂
∂ξi
(
A(ξ)Π
)
+
1
2
∑
ij
Bij
∂2Π
∂ξi∂ξj
(A.5)
where:
A(ξ) =
∑
j
Mijξj

Can a microscopic stochastic model explain the emergence of pain cycles in patients?13
with
M =
0
B
@
−α(1 − φ∗RT ) − βφ
∗
RT
− (γ + σ)φ∗H − δ1 − η1 c1/2[αφ∗T + β(1 − φ∗T − φ∗H)] −βφ∗RT − (γ + σ)φ
∗
T − η1
c1/2[α(1 − φ∗RT ) + βφ
∗
RT
] −c[αφ∗T + β(1 − φ∗T − φ∗H)] c1/2βφ∗RT
γφ∗H 0 γφ∗T − δ2
1
C
A
Notice that for c = 1 (see [14]), M reduces to the Jacobian matrix (10) which can be
directly calculated from the mean–field system. B is instead a symmetric matrix whose
elements reads:
B11 = αφ∗T (1− φ∗RT ) + βφ
∗
RT (1− φ
∗
T − φ∗H) + (γ + σ)φ∗Tφ∗H + δ1φ∗T + η1(1− φ∗T − φ∗H)
B12 = −c1/2[αφ∗T (1− φ∗RT ) + βφ
∗
RT (1− φ
∗
T − φ∗H)]
B13 = −γφ∗Tφ∗H
B22 = c[αφ∗T (1− φ∗RT ) + βφ
∗
RT (1− φ
∗
T − φ∗H)]
B23 = 0
B33 = γφ∗Tφ
∗
H + δ2φ
∗
H
References
[1] Victor M and Ropper A H. Principles Of Neurology. McGraw–Hill, 2001.
[2] Misulis K E and Fakhoury T. Spehlmann’s Evoked Potential Primer. Butterworth–Heinemann,
3rd edition, 2001.
[3] De Pascalis V and Cacace I. Pain perception, obstructive imagery and phase–ordered gamma
oscillations. Int. J. Psychophysiol., 56:157–169, 2005.
[4] Gross J, Schnitzler A, Timmermann L, and Ploner M. Gamma oscillations in human primary
somatosensory cortex reflect pain perception. PLoS Biol., 5:1168–1173, 2007.
[5] Rojas M J, Navas J A, and Rector D M. Evoked response potential markers for anesthetic and
behavioral states. Am. J. Physiol. Regul. Integr. Comp. Physiol., 291:R189–R196, 2006.
[6] Danhof M, de Jongh J, De Lange E C M, Della Pasqua O, Ploeger B A, and Voskuyl R
A. Mechanism–based pharmacokinetic–pharmacodynamic modeling: Biophase distribution,
receptor theory, and dynamical systems analysis. Annu. Rev. Pharmacol. Toxicol., 47:357–400,
2007.
[7] van Kampen N G. Stochastic preocesses in Physics and Chemistry. North Holland, Amsterdam,
1992.
[8] Grond S and Sablotzki A. Clinical pharmacology of tramadol. Clin. Pharmacokinet., 43:879–923,
2004.
[9] Gillen C, Haurand M, Kobelt D J, and Wnendt S. Affinity, potency and efficacy of tramadol and its
metabolites at the cloned human µ–opiod receptor. Naunyn–Schmiedeberg’s Arch. Pharmacol.,
362:116–121, 2000.
[10] Katzung B G. Basic & Clinical Pharmacology. McGraw-Hill, 9th edition, 2004.
[11] Lintz W, Barth H, Osterloh G, and Schmidt-Bothelt E. Bioavailability of enteral tramadol
formulations. 1st communication: capsules. Arzneimittelforschung, 36:1278–1283, 1986.
[12] McKane A J and Newman T J. Stochastic models in population biology and their deterministic
analogs. Phys. Rev. E, 70:041902, 2004.
[13] Gillespie D T. A general method for numerically simulating the stochastic time evolution of
coupled chemical reactions. J. Comp. Phys., 22:403–434, 1976.
[14] McKane A J and Newman T J. Predator-prey cycles from resonant amplification of demographic
stochasticity. Phys. Rev. Lett. , 94:218102, 2005.
[15] McKane A J, Nagy J D, Newman T J, and Stefanini M O. Amplified biochemical oscillations in
cellular systems. J. Stat. Phys., 128:165–191, 2007.

Can a microscopic stochastic model explain the emergence of pain cycles in patients?14
[16] Gardiner C W. Handbook of Stochastic Methods. Springer, second edition, 1985.
[17] Di Patti F and Fanelli D. In preparation.