## Abstract

In Lévy walks (LWs), particles move with a fixed speed along straight line segments and turn in new directions after random time intervals that are distributed according to a power law. Such LWs are thought to be an advantageous foraging and search strategy for organisms. While complex nervous systems are certainly capable of producing such behavior, it is not clear at present how single-cell organisms can generate the long-term correlated control signals required for a LW. Here, we construct a biochemical reaction system that generates long-time correlated concentration fluctuations of a signaling substance, with a tunable fractional exponent of the autocorrelation function. The network is based on well-known modules, and its basic function is highly robust with respect to the parameter settings.

## 1 Introduction

The migration paths of microorganisms can often be modeled as random walks. The simplest case of a (short-term) correlated random walk is the Ornstein-Uhlenbeck process, in which the autocorrelation function1 of the momentaneous velocity decays exponentially [19]. By contrast, in Lévy walks (LWs), particles move with a fixed speed, and turn abruptly after random time intervals that are distributed according to a power law (PL) [18], thereby leading to long-term correlated trajectories. Such LWs are thought to be an advantageous foraging and search strategy for the organisms [1,,,5, 7, 8, 11, 14, 16, 17, 20, 21]. However, it is not clear at present how single-cell organisms without complex nervous systems can generate the long-term correlated control signals required for a LW. Here, we show that such signals can arise simply by the combination of a Poisson point process with a positive feedback loop and a linear degradation mechanism, ingredients found frequently in even the simplest biochemical pathways [9].

We construct a chemical reaction system that generates random concentration peaks of a signaling substance, which could be used to trigger the different motor actions of a Lévy walker. We demonstrate that for a wide range of system parameters the intervals between trigger events are distributed according to a power law with fractional and tunable exponent γ, the hallmark of long-term correlations.

## 2 Concept

In what follows, we assume a well-stirred reaction system of unit volume, so that (continuous) concentrations and (integer) molecule numbers can be used interchangeably. We aim to design a reaction network that spontaneously creates a sequence of concentration pulses of a signaling substance B with a power law distribution p(τ) ∝ τ−γ of pulse durations. The autocorrelation of such a pulse train is long-time correlated with a fractional exponent in the range [−1, 0]. Concentration pulses of a given duration τ can be realized by quickly raising the concentration NB(t) of the signaling substance from NB(t1) ≈ 0 to a value NB(t2) ∝ τ that is proportional to the desired pulse duration, and thereafter decreasing NB(t) at a constant rate Rdeg, until it reaches zero again at time t3 (compare Figure 1).

Figure 1.

Sketch of the concentration b of signaling molecules B as a function of time t. Between times t1 and t2, concentration b is growing exponentially. After reaching a peak value NB(t2) it decays linearly to a value close to zero.

Figure 1.

Sketch of the concentration b of signaling molecules B as a function of time t. Between times t1 and t2, concentration b is growing exponentially. After reaching a peak value NB(t2) it decays linearly to a value close to zero.

In order to produce a power law distribution of initial concentrations NB(t2), we adopt a well-known mechanism [13]: If a positive random variable x is distributed as p(x) = ae−αx, then the function y(x) = beβx of this random variable is distributed according to
$py=abβy/b−1+αβ.$
(1)
The exponent of the power law distribution p(y) is, in general, a fractional number.

Using this scheme, we start from a small concentration NB(t1) of the signaling substance B at time t1 and let it grow exponentially. At time t2, the growth period of NB(t) is stopped by a Poisson point process. This way, the growth time t2t1 is exponentially distributed, and the combination with the exponential growth of NB(t) results in a power law distribution of the peak concentration NB(t2).

Chemically, the exponential growth of NB(t) is achieved via an autocatalytic reaction that has the property $ddtNBt∝NBt$. The Poisson point process is naturally realized by the diffusive arrival or disappearance of a diluted molecular species at the reaction site. The detailed chemical equations of our proposed reaction network are presented in the following section.

## 3 Detailed Reaction Network

### 3.1 Generation Reactions

The central chemical reaction of our proposed network is an autocatalytic reaction, in which molecules of chemical species B convert substrate molecules X into more B molecules, thus leading to an exponential growth of NB(t). In order to switch this autocatalytic process on and off, the reaction has additionally to depend on another coenzyme E:
$X+B+E→κ2B+E.$
(2)
If the supply of substrate X is large enough and the reaction is operating in the saturation regime, this results in an exponential growth of the number of B molecules, as can be seen from the following derivation.
The reaction's propensity function2 is given by
$π=κ⋅NE⋅NX⋅NB.$
(3)
Assuming that NXNB and NXNE for all times and therefore NX ≈ const, one can perform the substitution κ := κ ⋅ NX. Hence Equation 3 simplifies to
$π=κ⋅NE⋅NB.$
(4)
Since the net balance3 σB of molecule B is equal to 1, the temporal dynamics of NB is given by the differential equation
$N˙B=σB⋅π=κ⋅NE⋅NB$
(5)
$NBt∝eκ⋅NE⋅t.$
(6)

We have thus shown that NB is increasing exponentially during the growth period. If the duration t of this period could be randomly drawn from an exponential distribution, then at the end of the growth period the molecule number NB would have the desired power law distribution.

To generate exponentially distributed time periods t, we further assume that the coenzyme E can spontaneously decay into an inactive form:
$E→λ∅.$
(7)
Note that the decay rate of the coenzyme, affecting the parameter α in Equation 1, is proportional to the momentary molecule number NE. Therefore, starting with some initial amount, NE will decrease stepwise to zero. Each of these stepwise decays can be modeled as a Poisson point process with constant rate. It therefore follows that the time intervals t between the successive decay events are exponentially distributed:
$pt∝e−λ⋅NE⋅t.$
(8)
As the effective growth rate of B, corresponding to β in Equation 1, is also proportional to the number of coenzymes E, the final number NB(t2) of B molecules at the time t2 (when the last coenzyme E has disappeared) is the result of a series of distinct exponential growth periods. The final probability distribution NB(t2) is therefore not immediately obvious. However, since the power law exponent $−1+αβ$ in Equation 1 involves only the ratio of the coefficients α and β, the factor NE cancels out and the resulting distribution p(NB(t2)) is, indeed, the desired power law. The exponent of the power law can be expressed in terms of the rate constants λ and κ of the two chemical reactions as
$pNB∝NB−λ/κ+1.$
(9)

### 3.2 Decay Reactions

After the generation of power-law-distributed molecule numbers NB(t2), a further step is necessary to subsequently convert the molecule number into a duration τ with NB(t2) ∝ τ. This can be achieved by two coupled reactions: a synthesis reaction
$B+D→ν1F$
(10)
and a decay reaction
$F→ν2D+C.$
(11)
By design, the sum of the numbers of molecules of D and F is constant (ND+F = ND + NF = const). No new molecules of species D are generated, and all initially existing D molecules are either uncombined or part of the compound F. Taking into account that, due to exponential growth in the preceding reaction (2), most of the time NB is very large compared to the other molecule species involved, the propensity function
$π1=ν1⋅NB⋅ND$
(12)
of the reaction (10) is binary, that is, the quantity π1 is either equal to zero (in the absence of uncombined D) or very large. The latter is the case in the presence of at least one free D molecule. Thus a released D molecule is almost immediately recombined to another F. This means there are practically no free D molecules during the process (ND(t > 0) ≈ 0; compare Figure 2), and therefore
$NF=ND+F−ND≈ND+F=const.$
(13)
The propensity function
$π2=ν2⋅NF$
(14)
of the reaction (11) only depends on the reaction rate constant ν2 and the number of reactants, NF, which is constant as well. Hence, the temporal dynamics of NC is given by
$N˙C=π2=const$
(15)
$⇒N˙C=0.$
(16)
We have thus demonstrated that the generation rate of C is constant, so that NC is increasing linearly with time. This, in turn, means that NB is being degraded linearly with time (Figure 2). The degradation process comes to an end when no more B molecules are available and all compounds F have dissociated. This way, the total degradation time τ of the B molecules is proportional to the peak value of NB at the end of the exponential growth phase:
$τ∝NB.$
(17)
Considering the above-derived probability distribution p(NB), we find
$pτ∝τ−λ/κ+1.$
(18)
Figure 2.

Time course of molecule numbers NB, NC, ND, and NF during linear decay reactions. While NB decays linearly in time, NF grows linearly in time. There are practically no free D molecules during the process, and NC is nearly constant over time.

Figure 2.

Time course of molecule numbers NB, NC, ND, and NF during linear decay reactions. While NB decays linearly in time, NF grows linearly in time. There are practically no free D molecules during the process, and NC is nearly constant over time.

### 3.3 Switch Reactions

The four chemical reactions discussed above are not yet sufficient for proper function of repeated pulse generation.

On the one hand, the molecule species D disturbs the autocatalytic reaction (2): If some molecules B already decay by reacting with D during the generation phase, then NB does not grow exponentially.

On the other hand, D is absolutely indispensable once the decay reaction (7) has removed all enzymes E and hence has stopped the generation of B. In absence of D the previously described conversion process from NB to τ cannot start. Both problems can be tackled by a single modification of the reaction (7):
$E→D.$
(19)
In this case, the coenzyme E does not simply vanish but converts spontaneously into the species D. This reaction simultaneously stops the exponential generation and starts the linear decay of B. Thus, it acts as a switch between these two states.
The next problem arises when the two-stage decay process is finished, that is, when all B molecules have been converted to C (or BY) molecules. In this situation, the next cycle of generation and decay should begin. However, both essential ingredients to restart the autocatalytic reaction are missing, namely molecules of species E and B. To switch back from decay to production of B, we introduce a further reaction
$D+S→E$
(20)
where S is a coenzyme.
In order to start the autocatalysis, at least one B molecule in addition to the enzyme E is necessary. Hence we have to allow the production of B in another non-autocatalytic reaction:
$X+E→B+E.$
(21)

### 3.4 Input Reactions

In our reaction network, the species S acts as a coenzyme and represents the noisy input that is transformed into long-time correlated signals. Usually, coenzymes are continuously regenerated and their concentrations are maintained at a steady level inside the cell. We therefore may assume
$∅→S$
(22)
and
$S→∅.$
(23)
Here, the molecule species S can be regarded as a signaling substance.

### 3.5 Output Reactions

Finally, the reaction network is coupled to the motor unit of a model organism, which can produce either tumbling behavior (A1, during the short exponential growth phases of B) or running behavior (A2, during the linear decay phases). The coupling between the reaction network and the motor unit is formally described by the final two reactions:
$S→A1$
(24)
$C+S→A2$
(25)
where the specific rate constant κ9 is larger than κ8 (see below).

Since we have shown that the durations of the linear decay (running) phases are distributed according to a tunable power law, $pτ∝τ−αβ+1$ (compare Equation 1), this results in a Lévy walk of the organism.

### 3.6 Complete Reaction Network

Summing up, the complete reaction network consists of six molecule species (B, C, D, E, F, S) and the following ten chemical reactions:
$B+E→κ12B+EE→κ2B+EB+D→κ3FF→κ4C+DE→κ5DD+S→κ6E∅→κ7SS→κ8A1C+S→κ9A2S→κ10∅$

## 4 Simulation Results

The simulation of the chemical reactions has been implemented in C++ using the Gillespie algorithm [6]. The adjustable model parameters were both the initial numbers of molecules, Nn(t = 0), and the specific reaction rate constants κm. If not expressly described otherwise, the model parameters were set to the following default values:
$NB0=NE0=1NC0=ND0=NF0=NS0=0κ1=κ2=κ3=κ4=κ5=κ6=κ7=κ8=κ10=1.0κ9=5.0$

The time course of the molecule numbers NB, NC, ND, and NF during linear decay reactions is shown in Figure 2. While NB decays linearly in time, NF grows linearly in time. There are practically no free D molecules during the process, and NC is nearly constant over time.

Figure 3 shows the time course of molecule numbers NB, showing a series of peaks with a PL distribution of heights and widths. Each peak rises very fast and decays linearly with a constant rate, so that the ratio of height to width is the same for all peaks. To demonstrate the scale-invariance of the signal, we also show a zoom into the time course in Figure 4.

Figure 3.

Time course of molecule numbers NB, showing a series of peaks with a PL distribution of heights and widths. Each peak rises very fast (exponentially, not visible on this scale) and decays linearly with a constant rate. Hence, the ratio of height to width is constant for all peaks.

Figure 3.

Time course of molecule numbers NB, showing a series of peaks with a PL distribution of heights and widths. Each peak rises very fast (exponentially, not visible on this scale) and decays linearly with a constant rate. Hence, the ratio of height to width is constant for all peaks.

Figure 4.

Zoom into the first 1000 time units of the time course shown in Figure 3, demonstrating scale invariance.

Figure 4.

Zoom into the first 1000 time units of the time course shown in Figure 3, demonstrating scale invariance.

All above described results are very stable and highly robust to changes in rate constants. We generated 100 different reaction networks with modified rate constants by adding Gaussian-distributed random numbers with zero mean and standard deviation 0.2 to the reaction rates' default values to test the dependence of the essential functional modules on the parameter settings. In each case the resulting reaction networks yielded very similar results.

So far, we have shown analytically that the peak heights and peak widths of the fluctuating molecule number NB(t) are distributed according to a power law. If this is true, we would expect that the random time series s(t) = NB(t) has long-term correlations. To test this numerically, we have computed the autocorrelation function Css(τ) for different values of κ1 and κ5 (Figure 5). Indeed, we find that Css(τ) decays for long lag times as a PL with a fractional exponent that depends on the parameters κ1 and κ5.

Figure 5.

Autocorrelation function of the momentaneous molecule number s(t) = NB(t) of the signaling substance B for different values for κ1 and κ5. Blue: κ1 = 1.0 and κ5 = 1.0. Green: κ1 = 1.1 and κ5 = 1.0. Red: κ1 = 1.0 and κ5 = 1.2. All temporal correlations decay like a PL for large lag times, yet with different fractional exponents that can be tuned via the parameters κ1 and κ5.

Figure 5.

Autocorrelation function of the momentaneous molecule number s(t) = NB(t) of the signaling substance B for different values for κ1 and κ5. Blue: κ1 = 1.0 and κ5 = 1.0. Green: κ1 = 1.1 and κ5 = 1.0. Red: κ1 = 1.0 and κ5 = 1.2. All temporal correlations decay like a PL for large lag times, yet with different fractional exponents that can be tuned via the parameters κ1 and κ5.

## 5 Discussion

In the recent years, it has been increasingly acknowledged that biochemical reaction networks do not only work as deterministic machines, but that the stochastic nature of reactions, especially at low copy numbers of the molecules, can lead to useful functions for a cell [10, 15]. For example, it has been shown that covalent modification cycles can produce concentration fluctuations with Gaussian, exponential, or PL distributions and with tunable parameters [12].

In this work, we have demonstrated that also long-time correlated concentration fluctuations of a signaling substance can be generated by a system of relatively simple biochemical reactions. Coupling the concentration fluctuation to the motor units of an organism leads to scale-free random walks, which may be advantageous for efficient foraging and exploration. We therefore hypothesize that the essential functional modules of our reaction network (exponential concentration growth, linear concentration decay) may be identified in signaling pathways of single-cell organisms.

## Notes

1
The autocorrelation function of the time lag τ (more precisely, the autocorrelation coefficient) is defined as
$Cssτ=st+τ−y¯st−s¯tst−s¯2t$
where $s¯$ is the mean and 〈⋅〉t indicates averaging over time. Positive values of Css for a given τ mean that the values s(t) and s(t + τ) are similar with a high probability for all times t. Long-term correlated signals are marked by a power law decay of Css.
2
The propensity function of the mth reaction is defined as πm = κmγm(N), where κm is the specific rate constant and γm(N) is the number of distinct subsets of molecules that can form a reaction complex at time t, given by
$γmN=∏n∈NNnνnm.$
Nn is the number of molecules of the nth species at time t. The quantities νnm ≥ 0 are known as the stoichiometry coefficients of the reactants. These coefficients describe how many molecules of the nth species are consumed by the mth reaction.
3

The net balance σnm of the nth species in the mth reaction is defined as σnm = ν′nm − νnm, where νnm ≥ 0 and ν′nm ≥ 0 are the stoichiometry coefficients of the reactants and products, respectively.

## References

1
Bartumeus
,
F.
(
2007
).
Lévy processes in animal movement: An evolutionary hypothesis
.
Fractals
,
15
(
02
),
151
162
.
2
Bartumeus
,
F.
,
da Luz
,
M. E.
,
Viswanathan
,
G.
, &
Catalan
,
J.
(
2005
).
Animal search strategies: A quantitative random-walk analysis
.
Ecology
,
86
(
11
),
3078
3087
.
3
Benhamou
,
S.
(
2007
).
How many animals really do the Lévy walk?
Ecology
,
88
(
8
),
1962
1969
.
4
Codling
,
E. A.
,
Plank
,
M. J.
, &
Benhamou
,
S.
(
2008
).
Random walk models in biology
.
Journal of the Royal Society Interface
,
5
(
25
),
813
834
.
5
Edwards
,
A. M.
,
Phillips
,
R. A.
,
Watkins
,
N. W.
,
Freeman
,
M. P.
,
Murphy
,
E. J.
,
Afanasyev
,
V.
,
Buldyrev
,
S. V.
,
da Luz
,
M. G.
,
Raposo
,
E. P.
,
Stanley
,
H. E.
et al
, (
2007
).
Revisiting Lévy flight search patterns of wandering albatrosses, bumblebees and deer
.
Nature
,
449
(
7165
),
1044
1048
.
6
Gillespie
,
D. T.
(
1977
).
Exact stochastic simulation of coupled chemical reactions
.
The Journal of Physical Chemistry
,
81
(
25
),
2340
2361
.
7
Harris
,
T. H.
,
Banigan
,
E. J.
,
Christian
,
D. A.
,
,
C.
,
Wojno
,
E. D. T.
,
Norose
,
K.
,
Wilson
,
E. H.
,
John
,
B.
,
Weninger
,
W.
,
Luster
,
A. D.
et al
, (
2012
).
Generalized Lévy walks and the role of chemokines in migration of effector cd8+ t cells
.
Nature
,
486
(
7404
),
545
548
.
8
Humphries
,
N. E.
, &
Sims
,
D. W.
(
2014
).
Optimal foraging strategies: Lévy walks balance searching and patch exploitation under a very broad range of conditions
.
Journal of Theoretical Biology
,
358
,
179
193
.
9
Krauss
,
G.
(
2006
).
Biochemistry of signal transduction and regulation
.
Weinheim
:
Wiley-VCH
.
10
Levine
,
J.
,
Kueh
,
H. Y.
, &
Mirny
,
L.
(
2007
).
Intrinsic fluctuations, robustness, and tunability in signaling cycles
.
Biophysical Journal
,
92
(
12
),
4473
4481
.
11
Metzler
,
R.
,
Jeon
,
J. H.
,
Cherstvy
,
A. G.
, &
Barkai
,
E.
(
2014
).
Anomalous diffusion models and their properties: Non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking
.
Physical Chemistry Chemical Physics
,
16
(
44
),
24128
24164
.
12
Metzner
,
C.
,
Sajitz-Hermstein
,
M.
,
Schmidberger
,
M.
, &
Fabry
,
B.
(
2009
).
Noise and critical phenomena in biochemical signaling cycles at small molecule numbers
.
Physical Review E
,
80
,
021915
.
13
Newman
,
M. E.
(
2005
).
Power laws, Pareto distributions and Zipf's law
.
Contemporary Physics
,
46
(
5
),
323
351
.
14
Raichlen
,
D. A.
,
Wood
,
B. M.
,
Gordon
,
A. D.
,
Mabulla
,
A. Z.
,
Marlowe
,
F. W.
, &
Pontzer
,
H.
(
2014
).
Evidence of Lévy walk foraging patterns in human hunter–gatherers
.
Proceedings of the National Academy of Sciences of the U.S.A.
,
111
(
2
),
728
733
.
15
Rao
,
C. V.
,
Wolf
,
D. M.
, &
Arkin
,
A. P.
(
2002
).
Control, exploitation and tolerance of intracellular noise
.
Nature
,
420
(
6912
),
231
237
.
16
Reynolds
,
A. M.
(
2010
).
Bridging the gulf between correlated random walks and Lévy walks: Autocorrelation as a source of Lévy walk movement patterns
.
Journal of the Royal Society Interface
,
rsif20100292
.
17
Reynolds
,
A. M.
, &
Rhodes
,
C. J.
(
2009
).
The Lévy flight paradigm: Random search patterns and mechanisms
.
Ecology
,
90
(
4
),
877
887
.
18
Shlesinger
,
M. F.
,
Klafter
,
J.
, &
Wong
,
Y.
(
1982
).
Random walks with infinite spatial and temporal moments
.
Journal of Statistical Physics
,
27
(
3
),
499
512
.
19
Uhlenbeck
,
G. E.
, &
Ornstein
,
L. S.
(
1930
).
On the theory of the Brownian motion
.
Physical Review
,
36
(
5
),
823
.
20
Viswanathan
,
G.
,
Raposo
,
E.
, &
Da Luz
,
M.
(
2008
).
Lévy flights and superdiffusion in the context of biological encounters and random searches
.
Physics of Life Reviews
,
5
(
3
),
133
150
.
21
Viswanathan
,
G. M.
,
Buldyrev
,
S. V.
,
Havlin
,
S.
,
Da Luz
,
M.
,
Raposo
,
E.
, &
Stanley
,
H. E.
(
1999
).
Optimizing the success of random searches
.
Nature
,
401
(
6756
),
911
914
.

## Author notes

Experimental Otolaryngology, ENT-Hospital, Head and Neck Surgery, Friedrich-Alexander University Erlangen-Nuremberg (FAU), Germany. E-mail: holger.schulze@uk-erlangen.de

Department of Physics, Center for Medical Physics and Technology, Biophysics Group, Friedrich-Alexander University Erlangen-Nuremberg (FAU), Germany. E-mail: claus.metzner@biomed.uni-erlangen.de