Redundancy principle for optimal random search in

biology

Z. Schuss

∗

, K. Basnayake

†

, D. Holcman

†

Abstract

Chemical activation rate is traditionally determined by the diﬀusion ﬂux into

an absorbing ball, as computed by Smoluchowski in 1916. Thus the rate is set

by the mean ﬁrst passage time (MFPT) of a Brownian particle to a small target.

This paradigm is shifted in this manuscript to set the time scale of activation in

cellular biology to the mean time of the ﬁrst among many arrivals of particles at

the activation site. This rate is very diﬀerent from the MFPT and depends on

diﬀerent geometrical parameters. The shift calls for the reconsideration of physi-

cal modeling such as deterministic and stochastic chemical reactions based on the

traditional forward rate, especially for fast activation processes occurring in living

cells. Consequently, the biological activation time is not necessarily exponential.

The new paradigm clariﬁes the role of population redundancy in accelerating search

processes and in deﬁning cellular-activation time scales. This is the case, for ex-

ample, in cellular transduction or in the nonlinear dependence of fertilization rate

on the number of spermatozoa. We conclude that statistics of the extreme set the

new laws of biology, which can be very diﬀerent from the physical laws derived for

individuals.

1 Introduction

Why are specialized sensory cells so sensitive and what determines their eﬃciency? For

example, rod photoreceptors can detect a single photon in few tens of milliseconds [1],

olfactory cells sense few odorant molecules on a similar time scale, calcium ions can

induce calcium release in few milliseconds in neuronal synapses, which is a key process

in triggering synaptic plasticity that underlies learning and memory. Decades of research

have revealed the molecular processes underlying these cellular responses, that identify

molecules and their rates, but in most cases, the underlying physical scenario remains

unclear.

∗

Department of Applied Mathematics, Tel-Aviv University, Tel-Aviv 69978, Israel.

†

Computational Bioplogy and Applied Mathematics, IBENS, Ecole Normale Sup´erieure, Paris, France.

1

was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

The copyright holder for this preprint (whichthis version posted October 27, 2017. ; https://doi.org/10.1101/210443doi: bioRxiv preprint

Here, we discuss how the many molecules, which are obviously redundant in the tra-

ditional activation theory, deﬁne the in vivo time scale of chemical reactions. This redun-

dancy is particulary relevant when the site of activation is physically separated from the

initial position of the molecular messengers. The redundancy is often generated for the

purpose of resolving the time constraint of fast-activating molecular pathways. Activation

occurs by the ﬁrst particle to ﬁnd a small target and the time scale for this activation

uses very diﬀerent geometrical features (optimal paths) than the one described by the

traditional mass-action law, reaction-diﬀusion equations, or Markov-chain representation

of stochastic chemical reactions.

We also discuss the role of the fastest particle to arrive at a small target in the context

of fertilization, where the enormous and disproportionate number of spermatozoa, relative

to the single ovule, remains an enigma. Yet, when the number of sperms is reduced by

four, infertility ensues [2]. The analysis of extreme statistics can explain the waste of

resources in so many natural systems. So yes, numbers matter and wasting resources

serve an optimal purpose: selecting the most ﬁtted or the fastest.

More speciﬁcally, mass-action theory of chemical reactions between two reactants in

solution, A and B, is expressed as

A

k

1

k

−1

B,

where k

1

and k

−1

are the forward and backward reaction rates, respectively. The compu-

tation of the backward rate has long history that begins with Arrhenius’ law k

−1

= Ae

E/kT

,

where A is a constant, and E is activation energy, then Kramers’ rate, derived from the

molecular stochastic Langevin equation, which gives the prefactor A. For the past sixty

years, chemical physicists computed the activation energy E and clariﬁed the role of the

energy landscape, with extensions to applications in chemistry, signal processing (time

to loss of lock in phase trackers [3]), ﬁnance (time for binary option price to reach a

threshold), and many more [4].

In contrast, the forward rate k

1

represents the ﬂux of three-dimensional Brownian par-

ticles arriving at a small ball or radius a. Smoluchowski’s 1916 forward rate computation

reveals that

k

1

= 4πDca, (1)

where D is the diﬀusion coeﬃcient, when the concentration c is maintained constant

far away from the reaction site. When the window is in a smooth surface or inside a

hidden cusp, the forward rate k

1

is the reciprocal of the MFPT of a Brownian particle

to the window. The precise geometry of the activating small windows has been captured

by general asymptotics of the mean ﬁrst arrival time at high activation energy. This

mean time is, indeed, suﬃcient to characterize the rate, because the binding process is

2

was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

The copyright holder for this preprint (whichthis version posted October 27, 2017. ; https://doi.org/10.1101/210443doi: bioRxiv preprint

Poissonian and the rate is precisely the reciprocal of the MFPT. These computations are

summarized in the narrow escape theory [5–9]. For example, when the small absorbing

window represents binding at a surface, the forward rate is given by

k

−1

1

=

|Ω|

4aD

1 +

L(0) + N(0)

2π

a log a + o(a log a)

,

with |Ω| the volume of the domain of Brownian motion, a is the radius of the absorbing

window [5], and L(0) and N(0) are the principal mean curvatures of the surface at the

small absorbing window.

The forward rate k

1

has been used in almost all representations of chemical reactions:

it is the basis of the Gillespie algorithm for generating statistics of stochastic simulations

with rate k

1

. It has also been used in coarse-graining stochastic chemical-reactions into a

Markov chain. However, the rate k

1

is not used in simulations of Brownian trajectories.

In the latter case, the statistics of arrival times can be computed directly. Obviously,

diﬀusion theory computes k

1

from the mean arrival rate of a single particle. However,

in cell biology, biochemical processes are often activated by the ﬁrst (fastest) particle

that reaches a small binding target, so that the average time for a single particle do es

not necessarily represent the time scale of activation. Thus even the Gillespie algorithm

would give a rate, sampled with mean k

1

. But, as seen below, its statistics is diﬀerent

from the rate of the fastest arrival. This diﬀerence is the key to the determination of

the time scale of cellular activation that can be computed from full Brownian simulations

and/or asymptotics of the fastest particle. This is an important shift from the traditional

paradigm.

2 Biochemical reactions in cell biology

Chemical activation in cell biology starts with the binding of few molecules (Fig 1). The

signal is often ampliﬁed so that a molecular event is transformed into a cellular signal.

How fast is this activation? What deﬁnes its time scale? when there is a separation

between the site of the ﬁrst activation and that of ampliﬁcation. When particle move by

diﬀusion, is it suﬃcient that the ﬁrst particle arrives to a receptor to open it and lead to

an avalanche either through the entry of ions or the opening of the neighboring receptors.

Thus the time scale of activation is not given by the reciprocal of the forward rate, but

rather by the extreme statistics, that is, by the mean arrival time of the ﬁrst particle to

the activation site (the target).

The statistics of the ﬁrst particle to arrive to a target can be computed from the

statistics of a single particle when they are all independent and identically distributed

[10–13]. With N non-interacting i.i.d. Brownian trajectories (ions) in a bounded domain

Ω to a binding site, the shortest arrival time τ

1

is by deﬁnition

τ

1

= min(t

1

, . . . , t

N

), (2)

3

was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

The copyright holder for this preprint (whichthis version posted October 27, 2017. ; https://doi.org/10.1101/210443doi: bioRxiv preprint

RyR

First Ca

to arrive

2+

Activation of

an avalanche

mGluR

ER

IP3

to arrive

IP3

receptor

Ca

2+

Release

First four

Avalanche

arrive

AMPA/

First

by input currents

Neuro-

transmitter

release

Lost

receptor

(A)

(B)

(C)

Amplification

two to

molecules

Cell membrane

NMDA

Spine

Apparatus

Figure 1: (A) Calcium-induced-calcium-release in a dendritic spine. The ﬁrst Ryanodine Re-

ceptor (RyR) at the base of the spine apparatus that opens is triggered by the fastest calcium

ion. An avalanche of calcium release ensues by opening the neighboring receptors. This leads

to rapid ampliﬁcation at a much shorter time than the MFPT of the diﬀusing calcium ions.

(B) Activation of calcium release by IP3 receptors, which are calcium channels gated by IP3

molecules, which function as secondary messengers. When the ﬁrst IP3 molecules arrives at the

ﬁrst IP3R, its calcium release induces an avalanche due to the opening of subsequent IP3 recep-

tors. (C) In the post-synaptic terminal, the inﬂux of ions due to the opening of NMDA/AMPA

receptors is the ampliﬁcation process. The signal of the pre-synaptic signal transmitted by the

neurotransmitter molecules diﬀusing in the synaptic cleft and the time scale of the ampliﬁcation

are determined by the fastest molecules that arrive at the receptor targets and open them.

4

where t

i

are the i.i.d. arrival times of the N ions in the medium. The narrow escape

problem (NEP) is to ﬁnd the PDF and the MFPT of τ

1

. The complementary PDF of τ

1

,

Pr{τ

1

> t} = Pr

N

{t

1

> t}. (3)

Here Pr{t

1

> t} is the survival probability of a single particle prior to binding at the

target. This probability can be computed by solving the diﬀusion equation

∂p(x, t)

∂t

=D∆p(x, t) for x ∈ Ω, t > 0 (4)

p(x, 0) =p

0

(x) for x ∈ Ω

∂p(x, t)

∂n

=0 for x ∈ ∂Ω

r

p(x, t) =0 for x ∈ ∂Ω

a

,

where the boundary ∂Ω contains N

R

binding sites ∂Ω

i

⊂ ∂Ω (∂Ω

a

=

N

R

i=1

∂Ω

i

, ∂Ω

r

=

∂Ω − ∂Ω

a

). The single particle survival probability is

Pr{t

1

> t} =

Ω

p(x, t) dx, (5)

so that Pr{τ

1

= t} =

d

dt

Pr{τ

1

< t} = N(Pr{t

1

> t})

N−1

Pr{t

1

= t}, where Pr{t

1

= t} =

∂Ω

a

∂p(x,t)

∂n

dS

x

and Pr{t

1

= t} = N

R

∂Ω

1

∂p(x,t)

∂n

dS

x

. The pdf of the arrival time is

Pr{τ

1

= t} = NN

R

Ω

p(x, t)dx

N−1

∂Ω

1

∂p(x, t)

∂n

dS

x

, (6)

which gives the MFPT

¯τ

1

=

∞

0

Pr{τ

1

> t}dt =

∞

0

[Pr{t

1

> t}]

N

dt. (7)

The shortest ray from the source to the absorbing window δ

min

plays a key role, because

the fastest trajectory is as close as possible to that ray. The diﬀusion coeﬃcient is D,

the number of Brownian particles is N, s

2

= |x − A|, x the positions of injection, and

the center of the window is A. The asymptotic laws for the expected ﬁrst arrival time of

5