## 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 function^{1} 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 *N*_{B}(*t*) of the signaling substance from *N*_{B}(*t*_{1}) ≈ 0 to a value *N*_{B}(*t*_{2}) ∝ τ that is proportional to the desired pulse duration, and thereafter decreasing *N*_{B}(*t*) at a constant rate *R*_{deg}, until it reaches zero again at time *t*_{3} (compare Figure 1).

*N*

_{B}(

*t*

_{2}), 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

*p*(

*y*) is, in general, a fractional number.

Using this scheme, we start from a small concentration *N*_{B}(*t*_{1}) of the signaling substance B at time *t*_{1} and let it grow exponentially. At time *t*_{2}, the growth period of *N*_{B}(*t*) is stopped by a Poisson point process. This way, the growth time *t*_{2} − *t*_{1} is exponentially distributed, and the combination with the exponential growth of *N*_{B}(*t*) results in a power law distribution of the peak concentration *N*_{B}(*t*_{2}).

Chemically, the exponential growth of *N*_{B}(*t*) is achieved via an autocatalytic reaction that has the property $ddtNBt\u221dNBt$. 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

*B*convert substrate molecules

*X*into more

*B*molecules, thus leading to an exponential growth of

*N*

_{B}(

*t*). In order to switch this autocatalytic process on and off, the reaction has additionally to depend on another coenzyme

*E*:

*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.

^{2}is given by

*N*

_{X}≫

*N*

_{B}and

*N*

_{X}≫

*N*

_{E}for all times and therefore

*N*

_{X}≈ const, one can perform the substitution κ := κ ⋅

*N*

_{X}. Hence Equation 3 simplifies to

^{3}σ

_{B}of molecule

*B*is equal to 1, the temporal dynamics of

*N*

_{B}is given by the differential equation

We have thus shown that *N*_{B} 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 *N*_{B} would have the desired power law distribution.

*t*, we further assume that the coenzyme

*E*can spontaneously decay into an inactive form:

*N*

_{E}. Therefore, starting with some initial amount,

*N*

_{E}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:

*B*, corresponding to β in Equation 1, is also proportional to the number of coenzymes

*E*, the final number

*N*

_{B}(

*t*

_{2}) of

*B*molecules at the time

*t*

_{2}(when the last coenzyme

*E*has disappeared) is the result of a series of distinct exponential growth periods. The final probability distribution

*N*

_{B}(

*t*

_{2}) is therefore not immediately obvious. However, since the power law exponent $\u22121+\alpha \beta $ in Equation 1 involves only the ratio of the coefficients α and β, the factor

*N*

_{E}cancels out and the resulting distribution

*p*(

*N*

_{B}(

*t*

_{2})) 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

### 3.2 Decay Reactions

*N*

_{B}(

*t*

_{2}), a further step is necessary to subsequently convert the molecule number into a duration τ with

*N*

_{B}(

*t*

_{2}) ∝ τ. This can be achieved by two coupled reactions: a synthesis reaction

*D*and

*F*is constant (

*N*

_{D+F}=

*N*

_{D}+

*N*

_{F}= 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

*N*

_{B}is very large compared to the other molecule species involved, the propensity function

*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 (

*N*

_{D}(

*t*> 0) ≈ 0; compare Figure 2), and therefore

_{2}and the number of reactants,

*N*

_{F}, which is constant as well. Hence, the temporal dynamics of

*N*

_{C}is given by

*C*is constant, so that

*N*

_{C}is increasing linearly with time. This, in turn, means that

*N*

_{B}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

*N*

_{B}at the end of the exponential growth phase:

*p*(

*N*

_{B}), we find

### 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 *N*_{B} does not grow exponentially.

*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

*N*

_{B}to τ cannot start. Both problems can be tackled by a single modification of the reaction (7):

*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.

*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

*S*is a coenzyme.

*B*molecule in addition to the enzyme

*E*is necessary. Hence we have to allow the production of

*B*in another non-autocatalytic reaction:

### 3.4 Input Reactions

*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*can be regarded as a

*signaling*substance.

### 3.5 Output Reactions

*motor*unit of a model organism, which can produce either tumbling behavior (

*A*

_{1}, during the short exponential growth phases of

*B*) or running behavior (

*A*

_{2}, during the linear decay phases). The coupling between the reaction network and the motor unit is formally described by the final two reactions:

_{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\tau \u221d\tau \u2212\alpha \beta +1$ (compare Equation 1), this results in a Lévy walk of the organism.

### 3.6 Complete Reaction Network

*B*,

*C*,

*D*,

*E*,

*F*,

*S*) and the following ten chemical reactions:

## 4 Simulation Results

*N*

_{n}(

*t*= 0), and the specific reaction rate constants κ

_{m}. If not expressly described otherwise, the model parameters were set to the following default values:

The time course of the molecule numbers *N*_{B}, *N*_{C}, *N*_{D}, and *N*_{F} during linear decay reactions is shown in Figure 2. While *N*_{B} decays linearly in time, *N*_{F} grows linearly in time. There are practically no free *D* molecules during the process, and *N*_{C} is nearly constant over time.

Figure 3 shows the time course of molecule numbers *N*_{B}, 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.

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 *N*_{B}(*t*) are distributed according to a power law. If this is true, we would expect that the random time series *s*(*t*) = *N*_{B}(*t*) has long-term correlations. To test this numerically, we have computed the autocorrelation function *C*_{ss}(τ) for different values of κ_{1} and κ_{5} (Figure 5). Indeed, we find that *C*_{ss}(τ) decays for long lag times as a PL with a fractional exponent that depends on 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

_{t}indicates averaging over time. Positive values of

*C*

_{ss}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

*C*

_{ss}.

*m*th 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

*N*

_{n}is the number of molecules of the

*n*th 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

*n*th species are consumed by the

*m*th reaction.

The net balance σ_{nm} of the *n*th species in the *m*th reaction is defined as σ_{nm} = ν′_{nm} − ν_{nm}, where ν_{nm} ≥ 0 and ν′_{nm} ≥ 0 are the stoichiometry coefficients of the reactants and products, respectively.

## References

## 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