BPFD090623
BPFD090623
BPFD090623
H E S T U D I E S I T B E C A U S E H E D E L I G H T S I N I T, A N D H E D E L I G H T S
I N I T B E C A U S E I T I S B E A U T I F U L . I F N AT U R E W E R E N O T B E A U T I F U L ,
I T W O U L D N O T B E W O R T H K N O W I N G , A N D I F N AT U R E W E R E N O T
H E N R I P O I N C A R É , S C I E N C E A N D M E T H O D, 1 9 0 8
R AY M O N D E . G O L D S T E I N A N D E R I C L A U G A
BIOLOGICAL PHYSICS
A N D F L U I D DY N A M I C S
( A G R A D U AT E C O U R S E I N 2 4 L E C T U R E S )
D E PA R T M E N T O F A P P L I E D M AT H E M AT I C S A N D
THEORETICAL PHYSICS, UNIVERSITY OF CAMBRIDGE
Copyright © 2023 Raymond E. Goldstein and Eric Lauga
Preface 17
1 Overview (Lecture 1) 19
1.1 Length scales and methodologies 19
1.2 Scaling arguments 20
8 Bibliography 173
9 Index 183
List of Figures
4.18 The first three critical modes for the stretch-coil instability. 75
4.19 Experimental results on shear-induced buckling. 76
4.20 Setup for the follower force problem. 76
4.21 Components of complex eigenvalue for follower-force problem. 78
4.22 The leading eigenfunction of the follower-force problem. 78
4.23 Separable solution of follower-force problem. 79
Physics and Biology have been intertwined since the dawn of mod-
ern science, from Antony van Leeuwenhoek’s use in 1700 of ad-
vanced optics to reveal the hidden world of microscopic life,1 to 1
van Leeuwenhoek, A. IV. Part of a
Bonaventura Corti’s 1774 discovery of the persistent fluid motion letter from Mr. Antony van Leeuwen-
hoek, concerning the worms in sheeps
inside large eukaryotic cells,2 Robert Brown’s 1828 study of random livers, gnats, and animalcula in the ex-
motion at the microscale3 and Theodor Engelmann’s determination crements of frogs. Philos. Trans. R. Soc.,
22:509–18, 1700
in 1882 of the wavelength dependence of photosynthetic activity.4 2
B. Corti. Osservazioni Microscopische
Although at the time it might have been difficult to define the precise sulla Tremella e sulla Circolazione del Flu-
disciplines of each of these scientists—perhaps they were all ‘natural ido in Una Pianta Acquajuola. Appresso
Giuseppe Rocchi, Lucca, Italy, 1774
philosophers’—in hindsight we can see clearly the way in which their 3
Brown, R. XXVII. A brief account of
discoveries impacted both biology and physics. Despite this long his- microscopical observations made in the
tory of discoveries at the boundary between the two fields, and the months of June, July, and August, 1827,
on the particles contained in the pollen
innumerable fundamental contributions to both disciplines over the of plants; and on the general existence
long arc of time, the field of biological physics as a discipline within of active molecules in organic and in-
the research enterprise of physics has only risen to great prominence organic bodies. Philos. Mag., 4:161–73,
1828
since the postwar era, particularly since the mid 1980s. We are now 4
Engelmann, T.W. Ueber Sauerstof-
at the point that most academic physics departments have an iden- fausscheidung van Pfalnzenzellen im
Mikrospektrum. Pflüger, Arch., 27:
tifiable group in biological physics alongside those in high energy, 485–89, 1882
condensed matter, atomic and astrophysics.
As a subject, biological physics is in a much earlier stage of devel-
opment than the traditional ones found in physics or applied mathe-
matics curricula. This makes it an exciting field in which to work, but
also poses challenges for students and faculty, as there is no univer-
sally accepted canon to follow. Instead, each lecturer has had to cob-
ble together concepts and presentations from a number of disparate
sources to arrive at a syllabus. These lecture notes have been de-
veloped very much in that spirit. They derive from several master’s
level (Part III) courses taught by each of us, sometimes in collabora-
tion with Prof. Ulrich Keyser of the Cavendish Laboratory. We have
attempted to strike a balance between coverage, depth, and pedagogy
in order to address —in the 24 lectures of a Cambridge term—topics
in biological physics that cover length scales from the molecular to
the terrestrial. We presume no prior exposure to the subject, but
assume familiarity with the basic ideas of statistical physics, elec-
tromagnetism, nonlinear dynamics, and fluid mechanics. As this is
primarily a course for Part III students in the mathematical Tripos,
the presentation often includes digressions into the applied mathe-
matical aspects of the problems at hand.
We have organised the topics covered by these notes roughly in or-
18
der of increasing length scales, from the basic molecular forces that
govern the cohesion and interaction of molecules and supramolecu-
lar assembles, the elastic properties of filaments and membranes, the
role of thermal fluctuations in their conformation and interactions,
on to cellular motion, chemical kinetics and pattern formation. Our
approach is to give a sampling of the phenomena and their analy-
sis at each of the length scales of interest. Even with the goal of
brevity in mind, it is likely that covering all of these notes in 24 lec-
tures is a challenge. For further readings, we refer the student to
a number of excellent textbooks that have appeared in recent years
covering various aspects of biological physics5 and physical biology,6 5
P. Nelson. Biological Physics: Energy,
as well as more specialized works on elasticity.7 Much motivation for Information, Life. Chiliagon Science,
Philadelphia, PA, student edition, 2020
the physicist’s approach to biology comes from D’Arcy Wentworth 6
R. Phillips, J. Kondev, J. Theriot, H.G.
Thompson’s classic book,8 and we recommend highly the classic text Garcia, and N. Orme. Physical Biology of
on applied mathematics by Lin and Segel.9 Where possible, we have the Cell. Garland Science, Boca Raton,
FL, 2nd edition, 1998
included references to the original, often old literature, in the belief 7
B. Audoly and Y. Pomeau. Elasticity
that reading such scientific works is highly worthwhile.10 and Geometry: From hair curls to the non-
Associated with these notes is a set of Example Sheets. In many linear response of shells. Oxford Univer-
sity Press, Oxford, UK, 2010
places throughout the text we have deliberately left out certain cal- 8
D.W. Thompson. On Growth and Form.
culational details to encourage the reader to fill them in. At the Dover Publications, Mineola, N.Y., the
same time, we have taken pains to explain in a fair amount of detail complete revised edition, 1992
9
C.C. Lin and L.A. Segel. Mathematics
the ins and outs of various key calculations. While we emphasize Applied to Deterministic Problems in the
analytical calculations throughout the text, there are many contexts Natural Sciences. Macmillan Publishing
in which numerical work is an important adjunct to those calcula- Co., Inc., New York, 1974
10
Goldstein, R.E. Coffee stains, cell re-
tions, whether for something as simple as understanding the detailed ceptors, and time crystals: Lessons from
shape of a function, finding the roots of a transcendental equation, the old literature. Phys. Today, 71:32–38,
or to solve a nonlinear PDE to understand its global behaviour. To 2018
nm µm mm m km
hydrogen atom
matter. That is to say, gaining an intuitive feel for the typical sizes of
things, energy scales, and time scales is very important in building
a quantitative picture for a given problem. But what system of units
should we use? Many physicists and applied mathematicians have
been educated to use the MKS system, but its use can lead to very
awkward numbers; take for example the bending modulus k c of a
cell membrane, whose typical value might be 4 × 10−20 N·m, thus
having the units of energy. We shall see that the typical forces on the
cellular scale are picoNewtons (pN), and with typical sizes of molec-
ular constituents being nanometres (nm), the scale pN·nm is more
appropriate, for then k c ∼ 40 pN· nm. And, most importantly, with
the Boltzmann constant k B = 1.38 × 10−23 N·m/K and the absolute
temperature T = 300K, thermal energy kB T is ∼ 4 × 10−21 N·m or
k B T ≃ 4 pN · nm. (1.1)
Thus, concisely, the bending modulus is k c ∼ 10 k B T.
Table 1.1 collects the first set of physical quantities that will appear
in these notes. Other tables will follow as we develop the subject.
n2 h̄2 e2
E (r ) ∼ − . (1.2)
2mr2 r
The function E(r ) diverges to +∞ as r → 0 and approaches 0 from
below as r → ∞, having thus a single minimum. We find r by min-
imising E(r ), yielding the minimiser
n2 h̄2 me4
rn∗ = and En (rn∗ ) = − . (1.3)
me2 2n2 h̄2
r E
R= and E = , (1.4)
rn me4 /2n2 h̄2
1 2
E= 2
− , (1.5)
R R
which is minimised at R = 1, giving E = −1, as above. The ab-
sence of any explicit physical parameters from (1.5) indicates that a
"natural" system of units has been identified.
Let us now turn to an important biophysical question concerning
the role of inertia in the motion of microscopic organisms such as
bacteria.6 Suppose we have a spherical body of radius a and density 6
Purcell, E.M. Life at low Reynolds
ρ equal to that of the surrounding fluid, whose viscosity is µ. Let x (t) number. Am. J. Phys., 45:1–11, 1977
travel during that time? Without solving this explicitly we can sim-
ply assume there is a characteristic time τ and displacement ℓ (the
"coasting length") and balance terms, setting
4 3 ℓ ℓ
πa ρ 2 ∼ 6πµa (1.7)
3 τ τ
to obtain
2 a2
τ= and ℓ = v0 τ, (1.8)
9 ν
where ν = µ/ρ is the kinematic viscosity, with units of length2 /time.
From the structure of Eq. 1.6, we expect exponential decay of the daughter colonies
Setting these equal leads to the critical radius Figure 1.2: The green alga Volvox. It
consists of ∼ 103 biflagellated somatic
9ν
ac ∼ . (1.9) cells embedded in a transparent ex-
2v0 tracellular matrix (ECM), with a small
number of daughter colonies that grow
From this we see that even when the speed is ∼ 500 µm/s, which is up inside the spheroid.
among the highest speeds known for large multicellular algae such
as Volvox carteri shown in Fig. 1.2,7 the critical radius would still be 7
Goldstein, R.E. Green algae as model
∼ 1 cm, an order of magnitude larger than those organisms. Thus, organisms for biological fluid dynam-
ics. Annu. Rev. Fluid Mech., 47:343–375,
there is a vast range of organisms that swim in the low Reynolds 2015
number world. This subject will be covered in detail in Chapter 6.
We are now in a position to understand the relationship between
thermal energy and diffusion, as first done by Einstein8 and Suther- 8
Einstein, A. Investigations on the the-
land.9 The basic idea is to consider a suspension of microscopic ory of the Brownian movement. Ann. d.
Phys., 17:549–560, 1905
particles at concentration c( x ), with diffusion constant D, acted on 9
Sutherland, W. LXXV. A dynamical
by a force F ( x ) = −dϕ/dx, where ϕ is the potential energy. The flux theory of diffusion for non-electrolytes
J of particles is the sum of a Fick’s law contribution and an advective and the molecular mass of albumin.
Phil. Mag., 9:781–785, 1905
term with some speed u,
dc
J = −D + uc. (1.10)
dx
In the viscous regime, the speed u is determined by the force balance
ζu = −dϕ/dx, where ζ = 6πµa for spherical particles. To determine
D by the principles of equilibrium statistical physics we must indeed
be in equilibrium, with J = 0. Integrating (1.10), we obtain
c ∼ exp(−ϕ/Dζ ). (1.11)
kB T 0.2 µm2
D= ≃ . (1.12)
6πµa ã s
overview (lecture 1) 23
∂c
= D ∇2 c, (1.13)
∂t
we see by dimensional analysis that lengths ℓ and time t are related
by Dt ∼ ℓ2 or daughter colonies
ℓ2
t∼ . (1.14)
D
From the values of D derived above, we see that the time it takes
a small molecule to diffuse across a bacterium (ℓ ∼ 1 µm) is 10−3 s,
while that to cross a colony of Volvox (ℓ ∼ 300 µm) is about 103 s ∼
16 minutes. On the other hand, consider the aquatic plant Chara coral-
ECM
lina (Fig. 1.3). Its cells are some of the largest known to science;
they can reach ℓ ∼ 10 cm = 105 µm, giving a diffusive timescale of somatic cells 250 μm
p
≃ ρ + B2 ( T )ρ2 + B3 ( T )ρ3 + · · · . (2.2)
kB T
where the Bi are termed "virial coefficients". Intuitively, the quadratic TB T
term in ρ involves 2-body interactions, the cubic term captures 3-
body effects, and so on. Experimental measurements (Fig. 2.1) of the
second virial coefficient B2 show that it is negative at low tempera-
tures (indicating attraction between pairs) and becomes positive (re-
Figure 2.1: Temperature dependence of
pulsive) at high temperatures, crossing through zero at the so-called the second virial coefficient of a gas.
Boyle point TB .
One of the great early triumphs of statistical physics was to show
how B2 ( T ) arises from the underlying intermolecular potential. Its
derivation gives us an opportunity to review some basic features of
classical statistical physics that will be needed later. For the pur-
pose of this deviation, it is more convenient to work in the canonical
ensemble, although the more systematic approach uses the grand
26 biological physics and fluid dynamics
p2i
E= ∑ 2m + ∑ u(ri , rj ), (2.3)
i i< j
N2
Z
≃ 1+ d3 r f (r ) + · · · , (2.8)
2V
where, having done N − 1 integrals, there is one remaining factor
of V in the denominator. From the partition function Z we obtain
microscopic physics (lectures 2-4) 27
1
Z
Ea = Nρ d3 ru a (r ) , (2.10)
2
where again the factor of 1/2 avoids double-counting. We then de-
fine the parameter a as
1
Z
a=− d3 r u a (r ) , (2.11)
2
so Ea = − aNρ and the contribution to the pressure is −(∂Ea /∂V ) T,N =
− aρ2 . Thus, a first correction to the equation of state is
p ≃ ρk B T − aρ2 . (2.12)
Further, van der Waals realized that effect of the hard-core interac-
tions could be accounted for by subtracting from the total volume
V an amount proportional to the number of particles N, with an
effective excluded volume per particle b. Putting these two effects
together one has the revised equation of state
This has a form that is qualitatively like that seen in experiment (Fig.
2.1), with saturation at high temperatures to a constant, reflecting en-
tropic effects, and a divergence to −∞ at low temperatures as the at-
tractive part of the potential dominates. We conclude that the essence
28 biological physics and fluid dynamics
2πσ3 a
B2 ( T ) = − , (2.15)
3 kB T
where a is defined in (2.11), with the integration from σ to ∞. This
gives a precise meaning to the excluded volume parameter b; σ is
twice the hard-core radius, so b is four times the particle volume.
An important lesson from above is that interaction terms quadratic
in a density or concentration must be interpreted as free energies,
involving energy and entropy, rather than being purely energetic.
Moreover, we see how the integrated strength of the pairwise inter-
action, as appears in (2.11), arises naturally by mean field arguments.
They may also be understood as the zero-momentum values of the
Fourier transform of the potential, but for that transform to be well-
defined it is necessary to remove the singular repulsive contribution
to the potential, as we did in the decomposition of Fig. 2.2.
p21 1 p2 1
H0 = + mω02 x12 + 2 + mω02 x22 (2.16)
2m 2 2m 2
and the electrostatic interactions H1 ,
2 1 1 1 1
H1 = e + − − . (2.17)
R R − x1 + x2 R − x1 R + x2
microscopic physics (lectures 2-4) 29
2e2 x1 x2
H1 ≈ − , (2.18)
R3
x1 ± x2 x+ + x− x+ − x−
x± = √ ⇒ x1 = √ , x2 = √ . (2.19)
2 2 2
p21 1 2 2 p2 1 2 2
H≃ + mω+ x+ + 2 + mω− x− , (2.20)
2m 2 2m 2
where the new frequencies are
2 2e2
ω± = ω02 ± . (2.21)
mR3
Having diagonalised the Hamiltonian, the interaction potential u(r )
is just the change to the ground state energy of the system due to the
Coulomb interactions, namely the shift is zero-point energies,
1 1 1
u (r ) = h̄ω+ + h̄ω− − 2 · h̄ω0
2 2 2
2 2
(e /mω0 ) 2
1
≃ − h̄ω0 +··· , (2.22)
2 R6
where we have expanded the frequencies ω± to leading order. Cru-
cially, this is an attractive interaction. In grouping the factors in (2.22),
we have isolated the characteristic energy h̄ω0 that sets the overall
scale for the interaction. It would be typical of an internal excitation
energy from an s state to a p state, since the interaction is fundamen-
tally due to virtual transitions to states with dipole moments. Hence,
h̄ω0 is typically several electron volts, where 1 eV is 40 kB T.
Recognizing the R−6 dependence, we observe that e2 /mω02 must
be a characteristic volume (remember we are working in Gaussian
units!). To check this note that the simplest relationship between an
induced electric dipole moment d and an applied electric field E is
d = αE, (2.23)
1 1 eE0
H= mω02 z21 + mω02 z22 + · · · ; z1,2 = x1,2 + , (2.25)
2 2 mω02
1 h̄ω0 α2
u (r ) = − . (2.26)
2 r6
It follows that the typical scale of this energy within a condensed
phase of number density n is obtained by setting r −6 ∼ n2 , giving
8
Sculley, M.J., Duniec, J.T., Thorne,
termine the spacing observed in these structures has been analysed S.W., Chow, W.S. and Boardman, N.K.
along the lines we discuss below.8 The stacking of chloroplast thylakoids.
To calculate the interaction between slabs we start with the inter- Quantitative analysis of the balance of
forces between thylakoid membranes of
action of a single atom with a laterally infinite slab with atom num- chloroplasts, and the role of divalent
ber density n, as indicated in Fig. 2.5. Let us denote the interaction ions. Arch. Bioch. Biophys., 201:339–346,
1980
between two atoms as
C
u11 (r ) = − 6 . (2.28)
r
microscopic physics (lectures 2-4) 31
where we changed variables from z′ to z′′ = z − z′ in the second Figure 2.5: Geometry of an atom inter-
line (and dropped the prime). The power law r −6 of the interparticle acting with a laterally infinite slab of
thickness δ.
potential has become z−3 by integration over three dimensions.
Now we consider two slabs separated by a distance d (Fig. 2.6).
With the result 2.31, the slab-slab interaction uSS (d) is obtained by
integrating u1S (z) over the thickness of the second slab. If A is the
area of that slab, then there are Andz atoms in differential thickness δ
dz, and the interaction energy per unit area VA (z) ≡ uSS (z)/A is
Z d+δ d
VA (d) = dznu1S (z)dz (2.32)
d
δ
A 1 2 1
= − H − + , (2.33)
12π d2 (d + δ)2 (d + 2δ)2
Figure 2.6: Geometry of two slabs of
where have introduced the Hamaker constant9
thickness δ separated by face-to-face
distance d.
A H = π 2 n2 C, (2.34)
9
Hamaker, H.C. The London-van der
with units of energy. Based on our simple picture of the constant C Waals attraction between spherical par-
ticles. Physica, 4:1058–1072, 1937
in Eq. 2.26, we see that A H ∼ (1/2)π 2 h̄ω0 (αn)2 , so with h̄ω0 ∼ 3 eV
and αn ∼ 0.1 − 0.2 as discussed above, we estimate A H ∼ 5 − 25 kB T.
e2
E = eϕ(r ) = , (2.35)
r
where ϕ = e/r is the electrostatic potential. Essentially all of liv-
ing matter is infused with water, whose dielectric constant ϵ ∼ 80
32 biological physics and fluid dynamics
e2
E= . (2.36)
ϵr
While reduced by nearly two orders of magnitude, this interaction
still falls off as an inverse power of distance. A sense of the scale
of E is obtained by finding the distance at which it is comparable to
thermal energy. This Bjerrum length is
e2
λB = ∼ 0.7 nm. (2.37)
ϵk B T
Even pure water has small amounts of ionic species, since water
molecules are in equilibrium with protons (H+ ) and hydroxyl ions
(OH− ) according to the reaction scheme
H2 O ⇌ H+ + OH− . (2.38)
4π
∇2 ϕ = −
ϵ ∑ zs ec0s e−βzs eϕ . (2.41)
s
8πzec0
∇2 ϕ = sinh( βzeϕ). (2.42)
ϵ
This is a very nonlinear equation for which analytical solutions can
be found only in simplified geometries. Much of the important
physics of screened electrostatics can be seen in the weak field limit,
microscopic physics (lectures 2-4) 33
when βzeϕ ≪ 1, where we keep only the linearised r.h.s. of 2.42 and
obtain the Debye-Hückel equation
∇2 − κ 2 ϕ = 0, (2.43)
4πσ0 ϵκ
−n̂ · ∇ϕ = ⇒ σ0 = ϕ0 , (2.46)
z =0 ϵ 4π
4πσ0 −κz
ϕ(z) = e . (2.47)
ϵκ
Within Debye-Hückel theory the surface charge and surface poten-
tial are linearly related, but in the more general Poisson-Boltzmann
setting this relationship is nonlinear.
Observe that the DH equation ∇2 − κ 2 ϕ = 0 is the Euler-Lagrange
δG ∂ ∂(· · · ) ∂(· · · )
=− + , (2.49)
δϕ ∂x ∂ϕx ∂ϕ
The factors of 1/2 in Eqs. 2.50 and 2.51 arise from the linear relation
between surface potential and surface charge. In the more general
case, the free energy is expressed as a "charging integral" of the form
Z Z ϕ0
′ ′
− dS σ (ϕ )dϕ . (2.52)
S
G 2πσ02 F ϵκϕ02
= and =− , (2.53)
A ϵκ A 8π
where A is the total area of the surface.
cosh(κz)
ϕ = ϕ0 . (2.54)
cosh(κd/2)
microscopic physics (lectures 2-4) 35
Using this we find the charge density at the plate at z = d/2. Here,
−n̂ · ∇ = d/dz, thus
ϵκϕ0
σ (d/2) = tanh(κd/2) . (2.55)
4π
At the bottom plate −n̂ · ∇ = −d/dz, but with sinh(−d/2λ) =
− sinh(d/2λ) the charge density is the same. As the charge and
potential do not vary with position over these flat surfaces the sur-
face integration will just give a factor of the surface area A. The free
energy per unit area V (d) ≡ F (d)/A is
ϵκϕ02
V (d) = − tanh(κd/2). (2.56)
4π
In the limit d → ∞ we recover twice the single-surface result 2.53.
Subtracting the "self energy" of the surfaces, we obtain the repul-
sive interaction potential VR (d) = V (d) − V (∞),
ϵκϕ02
κd
VR (d) = 1 − tanh . (2.57)
4π 2
ϵκϕ02 −κd
VR (d) ≃ e , (2.58)
2π
which therefore has the same exponential decay as the electrostatic
potential itself. A second important point is that VR (d) is bounded
as d → 0; when the two charged surfaces are brought closer together,
their induced charges decrease fast enough to avoid a singularity.
Now we consider the case of two surfaces with fixed charge den-
sity σ0 . A simple calculation shows that the required potential is
4πσ0 cosh(κz)
ϕ( x ) = , (2.59)
κϵ sinh(κd/2)
4πσ02
κd
VR (d) = coth −1 . (2.60)
κϵ 2
ture that is a fusion of parts that are hydrophilic (having an affin- Figure 2.7: The chemical structure
of the phospholipid DPPC (dipalmi-
ity for water) and hydrophobic (repelling water). Figure 2.7 shows toylphosphatidylcholine).
the structure of the lung surfactant dipalmitoylphosphatidylcholine
(DPPC), which has two 16-carbon hydrophilic "tails" and a phos-
phatidylcholine head group. In diagrams of this type, the zig-zag
line represents the repeating group CH2 , terminated by the methyl
36 biological physics and fluid dynamics
group CH3 , with the hydrogen atom labels omitted for clarity. These
particular hydrocarbon chains are said to be "saturated" because ev-
ery carbon atom is 4-fold coordinated, bonded to two carbon atoms
and two hydrogen atoms (the maximum possible), unlike "unsatu-
rated" lipids with one or more C−C double bonds, thereby having
carbon atoms that have only a single bond to a hydrogen atom.
When placed in water, lipids self-organize into structures that iso-
late the oily hydrocarbon tails from contact with water molecules,
leaving the hydrophilic headgroups in contact with water. This "hy-
drophobic effect" plays a central role in cellular biology, dictating the 14
C. Tanford. The hydrophobic effect:
formation of micelles and biological mem-
structure of lipid assemblies, the folding of proteins, and localization branes. Wiley, New York, 1973
of channels in membranes.14 The simplest structure is a "micelle", a
water
spherical or cylindrical structure with the tails on the inside and the head groups
(but not always) the lipid molecules are free to diffuse laterally. A Figure 2.8: Structure of a lipid mem-
brane in water.
typical charged head group has a single elementary charge in an
area of ∼ 1 nm2 . If the Debye-Hückel screening length is ∼ 1 nm,
then 4πσ02 /κϵ ∼ 40 mN/m, which can be compared to the surface
tension of water (80 mN/m). Expressed another way, 40 mN/m is
40 pN·nm/nm2 , or about 10 kB T/nm2 .
Returning to the motivational problem of systems such as stacks
of thylakoid membranes, we now add the van der Waals interaction
2.33 between two membranes to the screened electrostatic repulsion
2.60. A convenient set of nondimensionalisations involves scaling
distances with the membrane thickness δ and the potential energy
with the Hamaker constant. Thus defining
d 12πδ2
x= , K = κδ, V= (VA + VR ), (2.61)
δ AH
we obtain
Kx 1 2 1
V ( x ) = Γ coth −1 − 2 + − , (2.62)
2 x ( x + 1)2 ( x + 2)2
where
24π 2 σ02 δ2
Γ= . (2.63)
ϵκA H
7 VA
-5
0 0.5 1 1.5 2
x
fixed enclosed volume of water less than the equivalent sphere (the
sphere having the same surface area as the membrane).
This functional is a surface integral of a quadratic form in the
curvatures,
1 1
Z
E = dS k c ( H − H0 )2 + kc K . (2.66) Figure 2.10: Electron micrograph of red
2 2 blood cells.
and based on the estimate 4πσ02 /ϵκ ∼ 10 kB T/nm2 above Eq. 2.61
and an assumed screening length of 1 nm, these moduli are on the
order of 10 kB T. Likewise, as the screening length is the only length
scale in the problem, we expect H0 ∼ κ −1 . The purpose of these
calculations is to make these estimates precise.
As we have already solved the problem for a single plane, we
next consider the cylinder. We wish to solve the modified Helmholtz
equation 2.43 inside and outside a cylinder of radius R, assuming
axisymmetry. The Debye-Hückel equation becomes
∂2
2 ∂ 2
r 2 + r − (κr ) ϕ = 0. (2.68)
∂r ∂r
4πσ I0 (κr )
ϕ (r ) = . (2.69)
ϵκ I1 (κR)
microscopic physics (lectures 2-4) 39
The (uniform) free energy per unit area of the inner surface is
1 1 4πσ02 I0 (κR)
σϕ = . (2.70)
2 S 2 ϵκ I1 (κR)
I0 (κR) 1 3
∼ 1+ + +··· . (2.72)
I1 (κR) 2κR 8(κR)2
ϕ( x, h( x )) = ϕ0 , (2.73)
1
ϕ( x, ϵh( x )) ≃ ϕ( x, 0) + εh( x )ϕz ( x, 0) + ε2 h( x )2 ϕzz ( x, 0) + · · · ,
2
(2.74)
where subscripts denotes partial differentiation. We expect that the
solution itself—in the bulk—also has an expansion in powers of ε,
1 (0) (1)
O ( ε2 ) : ϕ(2) ( x, 0) = − h2 ( x )ϕzz ( x, 0) − h( x )ϕz ( x, 0). (2.80)
2
A convenient way to solve these boundary-value problems is to
work in Fourier space. Let us define the Fourier transforms of the
potential and h with respect to x (keeping the z-dependence of ϕ),
Z Z
ϕ̂(m) (k, z) = dxeikx ϕ(m) ( x, z) and ĥ(q) = dxeiqx h( x ). (2.81)
and is solved by
ϕ̂(n) (q, z) = ϕ̂(n) (q, 0)e−κq z . (2.83)
To solve the sequence of boundary-value problems, we need the
Fourier transform of ϕ at z = 0; this can be determined directly from
the order-by-order boundary conditions. For example, at order ϵ1 we
have from (2.79)
ϕ̂(1) (q, 0) = κϕ0 ĥ(q). (2.84)
To complete the calculation, it is necessary to expand such quantities
as the surface normal vector
−h x êx + êz
n̂ = − p (2.85)
1 + h2x
e− βU
p( E) = , (3.2)
Z
and the expectation value of a quantity A (denoted by ⟨ A⟩) is
1
Z
⟨ A⟩ = dq N A(q N )e− βU . (3.3)
Z
When U is quadratic in the q N , as in the N = 1 case
1 2
U= kq , (3.4)
2
then R∞ 1 2 − βkq2 /2
−∞ dq 2 kq e
1 2 ∂ ln Z
kq = R∞ − βkq2 /2
=− . (3.5)
2
−∞ dqe
∂β
The final form on the r.h.s. of (3.5) tells us that ln Z is a generating
function. Now we change coordinates in Z to obtain
s
Z ∞ Z ∞
− βkq2 /2 2 2
Z= dqe = dxe− x , (3.6)
−∞ βk −∞
and thus
1
ln Z = − ln β + (terms independent of β). (3.7)
2
42 biological physics and fluid dynamics
To be clear, with the dual views of ensemble averages and time aver-
ages in statistical mechanics, this representation holds for any partic-
ular member of the ensemble of systems or at any particular moment
in time for a single string as it fluctuates; we view the amplitudes an
as the new fluctuating quantities and, in the usual way in statisti-
cal physics, we draw conclusions about the statistical properties of h
from those of the set { an }. By the orthogonality of the modes,
∞ nπ 2
γL
E=
4 ∑ L
a2n . (3.12)
n =1
2k B TL
⟨ a2n ⟩ = . (3.13)
γπ 2 n2
∞ ∞ mπx nπx
⟨h2 ( x )⟩ = ∑ ∑⟨ am an ⟩ sin sin
m n L L
∞
k B TL 2 sin2 (nπx/L)
=
γ π2 ∑ n2
n =1
kB T
= LF (ξ ), (3.14)
γ
where the final form defines the scaling function F shown in Fig. 3.2,
which depends on x and L only through the ratio ξ = x/L. The
linear scaling of the variance with system size L is reminiscent of
the scaling ℓ2 ∼ t for diffusion discussed in Chapter 1, and for the
analogous quantity in random walks, as we show later.
0.2
F (9)
0.1
0
0 0.5 1
9 = x=L
γL
E=
2 ∑ q2 |ĥ(q)|2 . (3.18)
q
The thermal average of the system average has a very compact form
in the case of the height function,
1 dq
Z
L ∑→ 2π
. (3.22)
q
kB T dq
Z
⟨ h2 ⟩ = . (3.23)
2πγ q2
In this case, the molecular cutoff scale a may be set to zero without
introducing any divergence in the average of interest, but, as men-
tioned earlier, the variance grows linearly with the large-scale cutoff
L and there is no sensible L → ∞ limit.
The average energy of the filament is proportional to the average
mean variance of the filament slope. A simple calculation shows
kB T
Z
γL 2
⟨ E⟩ = ⟨h ⟩ = L dq
2 x 4π
kB T L
= −1 . (3.25)
4 a
This result shows that the energy density ⟨ E⟩ diverges as the small-
1
Planck, M. Ueber das Gesetz der
scale cutoff is taken to zero. This "ultraviolet" divergence is exactly Energieverteilung im Normalspektrum.
what the quantum theory of blackbody radiation1 managed to tame Ann. Phys., 309:553–563, 1901
by introducing a discreteness to the modes at small scales.
If we generalize this calculation to d spatial dimensions and thus
a d − 1 dimensional surface, we find
d d −1 q k B T kB T q d −2
Z Z
⟨ h2 ⟩ ∼ ∼ dq (3.26a)
(2π )d−1 γq2 γ q2
k B T 1 q d −3 − q d −3 , d ̸ = 3
γ 3− d min max
∼ (3.26b)
∼ k B T ln ( L/a) , d = 3.
γ
and obtain
γA
E=
2 ∑(q2 + 2ℓ−c 2 )|ĥ(q)|2 , (3.29)
q
46 biological physics and fluid dynamics
This shows that the capillary length provides a cutoff on what would
otherwise be divergent fluctuation amplitudes as q → 0.
Introducing both small-scale and large-scale cutoffs 2π/a and 2π/L,
and assuming a circular geometry in Fourier space, the average vari-
ance of the displacement field is
" #
kB T qdq kB T 1 + 2 (π ℓc /a)2
Z
2
⟨h ⟩ = = ln . (3.31)
2πγ 2ℓ − 2
c +q
2 4πγ 1 + 2 (π ℓc /L)2
The thermodynamic limit L → ∞ is now possible at finite g, with
" 2 #
k B T π ℓ c
⟨ h2 ⟩ ∼ ln 1 + 2 . (3.32)
4πγ a
where ζ = 6πµa is the Stokes drag coefficient and ξ (t) is the ran-
dom force. For µm sized spheres and moderate laser power, k ∼ 10
fN/nm, so that when the particle is displacement is a fraction of the
particle size (1 µm = 103 nm), the trapping force will be in the pN
range. From (3.33) we see that there is a relaxation time
ζ
τ= . (3.34)
k
If a ∼ 1 µm and k is as above, τ ∼ 2 ms.
It is important to understand the meaning of the random force
ξ (t). First, note that the typical velocity of a water molecule is
√
found by equipartition to be 3k B T/m ∼ 600 m/s, and hence it takes
only about 3 × 10−13 s for a water molecule to traverse its own size
(0.2 nm), which in turn implies that the collision rate of molecules
with a test particle is perhaps 1020 s−1 . These time scales are vastly
shorter than anything we measure when tracking the microsphere’s
motion, and it is thus natural to exploit this separation of time scales
in the analysis. A second point regards what it means to "solve" a
stochastic ODE such as (5.41). As we shall see below, for any par-
ticular realisation of the random noise ξ (t) the Langevin equation is
simply a linear ODE with a time-dependent force. Yet, in the typical
experiment we perform averages over many stochastic trajectories
of the particle to arrive at the average behaviour. This process cor-
responds to averaging over the realisations of ξ (t). Denoting such
averages with angular brackets, we assume ⟨ξ ⟩ = 0 by symmetry.
To see these issues in action, first simplify (3.33) to
ẋ + τ −1 x = η (t), (3.35)
If we average this result over realisations of the noise, then the deter-
ministic term is unchanged, while the integral term vanishes by the
assumed zero value of the mean of the noise, leaving
Γτ
⟨ x ( t )2 ⟩ = . (3.41)
2
Here we appeal to equipartition, since the linear force law (3.33)
arises from an energy kx2 /2, which in turn implies ⟨ x2 ⟩ = k B T, so
2k B T
Γ= . (3.42)
ζ
kB T
D= (3.44)
ζ
0.5
0
0 1 2 3 4 5
f
rn +1 = rn + ζ n , (3.56)
N
R ≡ r N − r0 = ∑ ζn (3.57)
n =1
N
⟨R⟩ = ∑ ⟨ζ n ⟩ = 0 (3.58)
n =1
1
P({rk }) = G ({rk }) with G = e− βU ({rk }) . (3.60)
Z
The energy is typically a sum of near-neighbor interactions and a
contribution from an external potential,
N
U ({rk }) = ∑ Uj (rj−1 , rj ) + W ({rk }). (3.61)
j =1
d3 k ik·P
Z
δ (P) = e , (3.64)
(2π )3
to express G for a system of identical segments as
d3 k −ik·R
Z
G= e K (k; N ), (3.65)
(2π )3
where the so-called "characteristic function" is
Z N
K (k; N ) = dR p(R)eik·R . (3.66)
3R2
Fentropic ( R) = −k B T ln G( R; N ) = k B T . (3.71)
2Nb2
This has the form of a (Hookean) entropic spring, with a minimum at
R = 0 reflecting the fact that stretching the chain lowers its entropy.
In a root-mean-square sense, the typical end-to-end distance obeys
the scaling form
R ∼ Nν, (3.72)
with ν = 1/2 as expected for an ideal random walk. Let us now
consider how excluded-volume interactions change this free energy,
fluctuations & brownian motion (lectures 5-7) 53
N vN 2
vk B T · N · ∼ k B T . (3.74)
R3 R3
The total free energy in this Flory theory is thus
3R2 vN 2
F ( R) = k B T + . (3.75)
2N ℓ2 R3
R∗ ∼ N 3/5 . (3.76)
3
ν= . (3.78)
d+2
These are compared to the known values in Table 3.1. The case of
d = 1 is trivial, as such a self-avoiding chain must be fully stretched
out, while for d > 4 the probability of a chain contacting itself is so
low as to leave the Gaussian scaling unchanged. Apart from the case
d > 4 we see that Flory theory is remarkably accurate.
Living systems are replete with filaments, from the polymers such
as actin and microtubules that comprise the cytoskeleton to cilia and
flagella that power the motility of microorganisms. At larger scales
we find slender, undulating organisms such as worms and snakes. In
this chapter we discuss the basic differential geometry of filaments
in two and three dimensions, their elasticity and principles by which
their motion can be understood. We also discuss some key insta- 1
Powers, T.R. Dynamics of filaments
bilities that are relevant in biological phenomena. Aspects of the and membranes in a viscous fluid. Rev.
analyses presented here find application in the theory of pattern for- Mod. Phys., 82:1607–1631, 2010
mation, discussed in Chapter 7. There are many excellent reviews on 2
S. Lipschutz. Differential Geometry.
the subject,1 and concise compendia on differential geometry.2 Schaum’s Outline Series. McGraw-Hill
Education, New York, 1969
where ϵij is the Levi-Civita symbol (ϵ12 = +1, ϵ21 = −1, ϵ11 = ϵ22 =
0). The pair (t̂, n̂) is a local coordinate system (Fig. 4.2) that rotates
56 biological physics and fluid dynamics
where T and V are the kinetic and potential energies. From the
R
action S = dtL, the minimum action principle δS = 0, yields the
Newtonian equation of motion
d ∂L ∂L
− = 0. (4.14)
dt ∂q̇ ∂q
∂V
mq̈ = − . (4.15)
∂q
In the limit of vanishing inertia (Re = 0), simply dropping the kinetic
energy term leaves us with no dynamics, because dissipative forces
are not captured in either T or V . Instead, we must introduce a gener-
alised force associated with dissipation. This quantity is related to the
rate at which viscous forces do work. For example, if, as in the case
of a sphere moving through a viscous fluid, the viscous force is ζ q̇,
where ζ is a drag coefficient, the rate of dissipation (force×velocity)
is ζ q̇2 . Thus, we introduce the Rayleigh dissipation function
1 2
R= ζ q̇ , (4.16)
2
proportional to the dissipation rate, and obtain the new variational
principle in which the derivative of R with respect to velocity is a
generalised force,
d ∂L ∂L ∂R
− =− . (4.17)
dt ∂q̇ ∂q ∂q̇
In the overdamped limit, we then have the Aristotelian law
∂V
ζ q̇ = − . (4.18)
∂q
58 biological physics and fluid dynamics
1 δE
ζrt = − √ , (4.20)
g δr
where I = t̂t̂ + n̂n̂ + b̂b̂ is the unit dyad. If we use the asymptotic
result ζ ⊥ = 2ζ ∥ , then the projection operator on the l.h.s. of 4.23 is
proportional to I − t̂t̂/2. The inverse of this operator is
−1
1
I − t̂t̂ = I + t̂t̂, (4.24)
2
and therefore we can write the equation of motion 4.23 as
1 1 δE
rt = − I + t̂t̂ · √ . (4.25)
ζ⊥ g δr
m = ζ r ω, (4.26)
4.2.2 Elasticity
In this section we use the general formulation developed above to
study filaments governed by the elastic energy7 7
L.D. Landau and E.M. Lifshitz. Theory
of Elasticity. Butterworth-Heinemann,
Z L
1 Oxford, UK, 3rd edition, 1986
E= A ds κ 2 , (4.27)
2 0
flexible", and are found throughout the cell. Table 4.1 summarizes
the three most important examples. We shall not delve into the de- extension
tailed derivation of 4.27 here, but note that resistance to bending
fundamentally arises from the resistance to stretching. As shown in
Fig. 4.3, the outer edge of a bent filament is under extension and the compression
inner edge is compressed. The resistance to stretching or compres-
Figure 4.3: A bent filament viewed in
sion is quantified by the Young’s modulus Y as the proportionality longitudinal cross section. Red dashed
between the force per unit area applied to an object and its fractional line indicates neutral plane.
change in length,
force ∆L
=Y . (4.28)
area L
Y clearly has units of energy/volume, and by comparison to the
modulus we deduce that A ∼ Y · length4 . As the only length scale of
relevance for a filament is its radius a, we deduce that
A ∼ Ya4 . (4.29)
√
The quantities rα and rαα show up in the metric factor g and the
elastic term, which can be written in a covariant manner as
2
1 rα
κ 2 = √ ∂α √ . (4.31)
g g
δE ∂f ∂ ∂f ∂2 ∂ f
= − + 2 +··· , (4.32)
δr ∂r ∂α ∂rα ∂α ∂rαα
where each successive term arises from one additional integration by
parts in the usual manner.
For the remainder of our discussion we will confine ourselves to
curves in the plane. A straightforward (if tedious) calculation of this
to 4.27 for such filaments yields the result
1 δE
1 3
−√ = A κss + κ n̂. (4.33)
g δr 2
Note that this force is purely normal to the filament. As the curva-
ture involves two derivatives of position, κss involves four; this is a
hallmark of elasticity that is distinct from the effects of tension (as
in a violin string), which produces a restoring force proportional to
the curvature itself. One might imagine that with an energy that
is (curvature)2 the fourth derivative would be the only term in the
force, but this is not the case. The nonlinear term κ 3 appears even
though this theory is based fundamentally on linear elasticity the-
ory. It is a geometric nonlinearity that arises from the fact that the
coordinate r and arclength s are not independent. And it is intuitive
physically, if we think of a closed circular loop, where the curvature
is constant (κss = 0) and yet there is uniform outward force which
would only be relaxed at vanishing curvature (infinite radius).
A small technical aside is needed here. It is tempting (and fre-
quently done in the literature) to adopt a simplified picture in which
the energy (4.27) as a functional of the quantity r(s), ignoring the fact
that s and r are not independent. In this calculation, we simply write
1
Z
E [r] = A ds r2ss , (4.34)
2
and deduce the force
∂ 1 2
−∂ss Ar = − Arssss . (4.35)
∂rss 2 ss
1 δE
3 2
−√ = − Arssss − ∂s Aκ t̂ . (4.37)
g δr 2
biological filaments (lectures 8-12) 61
We shall now see that the difference between the two calculations
represents a tension-like contribution.
If we were to use E directly as our energy functional, the dynamics
will not necessarily preserve the local arclength along the filament.
To enforce local inextensibility we introduce a Lagrange multiplier
function Λ(α) that represents the internal tension,8 and define 8
Goldstein, R.E. and Langer, S.A. Non-
linear dynamics of stiff polymers. Phys.
Z 1
√ Rev. Lett., 75:1094–1097, 1995
F =E− dα gΛ(α). (4.38)
0
rt = U n̂ + W t̂, (4.40)
We now see that the difference between the “proper" (4.33) and sim-
plified (4.36) forces is a redefinition of the tension.
The general dynamics 4.40 is an intrinsic equation of motion, and it
is of interest to the time evolution of geometric quantities such as the
tangent angle θ and the curvature κ. As these are scalar quantities
they are often easier to study (both analytically and numerically),
particularly when a filament is free in space and its absolute position
62 biological physics and fluid dynamics
θt = −Us + θs W. (4.44)
A bit more work shows that the curvature obeys the PDE
κt = − ∂ss + κ 2 U + κs W . (4.45)
AW4x = k4 W, (4.50)
-0.5
-1
0 : 2: 3:
kL
-1
-2
0 0.5 1
x=L
h( x ) = ∑ a n W ( n ) ( x ), (4.54)
n
where in the last relation λn is the nth eigenvalue, and we have used
the orthonormality the Ws to simplify the result. Once again we
have found an energy expressed as a sum of independent quadratic
contributions, and thus from equipartition we obtain
kB T k T 16L4
⟨ a2n ⟩ = ≃ B , (4.56)
Aλn A (2n + 1)4 π 4
where in the last relation we have used the approximate form of the
roots of the transcendental equation (4.53). The variance is
kB T W ( n ) ( x )2
⟨h2 ( x )⟩ =
A ∑ λn
. (4.57)
n
imposed flow fields was done some time ago, using a microfluidic
setup depicted in Fig. 4.6.10 Single actin filaments can be held at 10
Kantsler, V. and Goldstein, R.E. Fluc-
the crossing point of two perpendicular channels meeting at a four- tuations, dynamics, and the stretch-coil
transition of single actin filaments in ex-
fold junction by means of the pressure difference ∆p between the tensional flows. Phys. Rev. Lett., 108:
inlets and outlets. This also allows the flow to be changed from com- 038103, 2012
pressional to extensional for a filament oriented along either of the
orthogonal directions. The vector field show illustrates the hyper-
bolic flow profile that can be achieved. At right is an image of a
fluorescently labelled actin filament and the tracing of its centreline.
The variance of the filament deviation h( x ) from its straight con-
figuration, normalised by its value at the endpoint, is shown in Fig.
4.7 for these experiments with actin, compared with the contribution
biological filaments (lectures 8-12) 65
11
Fygenson, D.K., Marko, J.F. and
Libchaber, A. Mechanics of micro-
tubule-based membrane extension.
Inside of cells, microtubules often interact with cell membranes Phys. Rev. Lett., 79:4497–4500, 1997
and deform them. A series of experiments in the 1990s11 provided
some of the nicest illustrations of the forces at work in these con-
texts. Figure 4.8 shows an experimental realization of this instability
involving microtubules grown by inducing polymerisation of tubu-
lin monomers in solution inside of lipid vesicles. As they grew and
the deformed membrane exerts larger and larger thrusting forces on y
their ends, the microtubules buckle.
Our analysis of Euler buckling will be done using the tangent- θ s
angle representation of the curve to illustrate an alternate approach h(x)
to the Monge representation. Using the geometry of Fig. 4.9, the F x
coordinates x (s) and y(s) of a point at arclength s are related to Figure 4.9: Geometry of Euler buckling.
the tangent angle through dx/ds = cos θ and dy/ds = sin θ. If we
66 biological physics and fluid dynamics
take the left end of the filament to be at the origin (x (0) = 0), then
the difference between the fully extended length L and the projected
distance between the endpoints is x ( L) and x (0) is
Z L Z L
dx
L − x ( L) = ds 1 − = ds (1 − cos θ (s)) . (4.58)
0 ds 0
The energy functional for the filament includes the bending energy
and the work done by the force F in displacing the filament ends,
Z L
1 2
E[θ ] = ds Aθs − F (1 − cos θ ) . (4.59)
0 2
FL2
θαα + f sin θ = 0, with f = . (4.61)
A
Let us consider clamped boundary conditions, θ (0) = θ (1) = 0.
From everyday experience we expect the buckled shape beyond some
threshold force to have a single maximum in the middle of the fila-
ment, with θ = 0 there as well. This suggests a trial function
θαα ≃ − f θ. (4.63)
f = f c (1 + ϵ ) , (4.65)
2 0.5
(a) (b) Figure 4.11: Euler buckling. (a)
Numerically obtained solutions θ (s)
1
for ϵ = 0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.7,
with clamped boundary condi-
3(s)
y=L
0
tions. (b) The corresponding
filament shapes. Obtained with
-1
Figure_411_and_412_Eulerbuckle.m.
-2 0
0 0.5 1 0 0.5 1
s=L x=L
Figure_411_and_412_biharmonic.m.
1
0
0 0.5 1
0
La2
Z L
1
∆x ≃ ds θ 2 ≃ +··· (4.68)
0 2 4
68 biological physics and fluid dynamics
√
As a = 8ϵ near the bifurcation, we obtain
Fc
F − Fc = ∆x, (4.69)
2L
and thus beyond the buckling the effective spring constant is Fc /2L.
For Stokes Problems I and II (SI & SII), we start from the Navier-
Stokes equations for an incompressible fluid with density ρ, viscosity
µ, pressure p and velocity field u,
and recognize for wall-driven flow that the nonlinear term vanishes
by symmetry, and we can ignore the pressure gradient. If the wall
moves along êy , then the velocity will be of the form u = u( x, t)êy . It
follows that the function u obeys the diffusion equation
ut = νuyy , (4.71)
ζ ⊥ ht = − Ah xxxx . (4.77)
We focus on the case of EHDII and, as with SII, consider the situa-
tion post-transients in which the left-hand side is forced as h(0, t) =
h0 cos(ωt), with zero torque (h xx (0, t) = 0), and the distant end is
free, with h xx ( L) = h xxx ( L) = 0. By analogy with SII, dimensional
analysis shows there is an elastohydrodynamic length
1/4
A
ℓ(ω ) = , (4.78)
ζ⊥ω
70 biological physics and fluid dynamics
λ1 = C8 − iS8 , λ2 = S8 + iC8 ,
λ3 = −C8 − iS8 , λ4 = −S8 − iC8 , (4.81)
-1
0 2 4 6
@
-1
1
2
0
-1
1
L=1
0
-1
0 1 2 3 4 5 6 7
@
the "scallop theorem". Since the Stokes equation 0 = −∇ p + µ∇2 u number. Am. J. Phys., 45:1–11, 1977
movement to the right at an equal and opposite speed. But the rod’s
motion is indistinguishable in the two cases, so the only allowable
speed is zero. This is referred to as the “scallop theorem" as a means
of highlighting that an object with a single degree of freedom exe-
cuting cyclic motion will exhibit no net propulsion. As the rod starts
to flex during the oscillations, time reversal invariance is broken, so
the two movies are not equivalent, and net propulsion can occur.
EHDII provides a particularly simple way to understand this ef-
fect, by means of the time-averaged x-component of the force asso-
ciated with the filament oscillations. Here, we project the left hand
side of the general force law (4.23) on to the x axis, yielding for our
problem of a curve in the plane17 17
Camalet, S. and Jülicher, F. and Prost,
J. Self-organized beating and swim-
F = −êx · ζ ∥ t̂t̂ + ζ ⊥ n̂n̂ · rt ≃ ζ ⊥ − ζ ∥ h x ht + · · · . (4.84) ming of internally driven filaments.
Phys. Rev. Lett., 82:1590–1593, 1999
where we have used the leading order tangent and normal vectors
and written rt ≃ ht êy . Integrating this over the filament yields the
total propulsive force generated. Using (4.77), we have
Z L Z L
A
dxh x ht = − dxh x h4x . (4.86)
0 ζ⊥ 0
1
⟨ Fp ⟩ = ζ ⊥ − ζ ∥ h20 ωY (L), (4.90)
2
where L = L/ℓ(ω ) and
0.3
Z
0.2
0.1
0
0 5 10
L
with F (χ) = F R (χ) + iF I (χ). The function Y (L) is shown as the blue
curve in Fig. 4.17. The propulsive force indeed tends to zero as the L
decreases and the shapes are time-reversal-invariant, consistent with
the scallop theorem. As L increases, Y displays a maximum followed
by attenuated oscillations as it asymptotes to a finite value as L → ∞.
This saturation indicates that beyond L ∼ 8 no more of the filament
is set in motion and no further propulsive force is gained.
In the regime where the propulsive force has saturated any ad-
ditional length added to the filament serves only to increase the
drag for motion along the x axis, and one therefore expects that
the propulsion speed of the filament would peak at an intermediate
value of L. As we have assumed in EHDII that the filament is nearly
straight (|h x | ≪ 1, it is consistent to assume that the average lateral
propulsive speed ⟨v x ⟩ of the undulating filament can be deduced by
the simple balance of force ζ ∥ L⟨v x ⟩ = ⟨ Fp ⟩, or
⟨ Fp ⟩ h2 ω
⟨v x ⟩ = = 0 V (L), (4.92)
ζ∥ L 2ℓ(ω )
where V (L) = Y (L)/L. As seen in Fig. 4.16, the peak in ⟨ Fp ⟩ implies
a peak in speed, and ⟨v x ⟩ vanishes as L → ∞ as anticipated.
The average power P∥ associated with the lateral motion is P∥ =
v x ⟨ Fp ⟩ ∼ ⟨ Fp ⟩2 . It can be compared to the power
Z L
P⊥ = ζ ⊥ ⟨ dxh2t ⟩ (4.93)
0
generated by transverse motions. Together, these define the "effi-
ciency" ε = P∥ /P⊥ of the propulsive mechanism, which has the form
2
h0
ε∼ Z (L), (4.94)
ℓ(ω )
where the scaling function is
Y 2 (L)
Z (L) = RL . (4.95)
L 0 dχ| F (χ)|2
Figure 4.17 shows that the efficiency also displays the peaked be-
haviour of the other scaling functions, suggesting an optimum length.
74 biological physics and fluid dynamics
σx = −ζ ∥ γ̇x, (4.96)
1 L2
σ( x) = − x2 ζ ∥ γ̇, (4.97)
2 4
a relation noted much earlier by Batchelor.19 In the presence of this 19
Batchelor, G.K. Slender-body theory
nonuniform tension, the filament energy functional is for particles of arbitrary cross-section in
Stokes flow. J. Fluid Mech., 44:419–440,
Z L/2 1970
1
E= dx Ah2xx + σ( x )h2x . (4.98)
2 − L/2
µγ̇L4
Σ= √ , (4.102)
π 3 ln( L/a e)
biological filaments (lectures 8-12) 75
and the logarithm in the denominator arises from the RFT drag coef-
ficient. Note that the tension term on the r.h.s. of (4.101) is not a total
derivative; the linearised equation of motion is not variational. This
arises from the anisotropic drag and the fact that the background
flow that enters the drag force is the source of the tension.
The determination of the threshold for the stretch-coil instability
is a matter of a numerical solution of the the eigenvalue problem that
comes from a separable solution H (χ, τ ) = exp(ωτ ) F (χ), namely
2
π 2
F4χ − Σ − χ Fχχ − 4χFχ = 4Σ (ω − 1) F. (4.103)
4
Our general expectation is that there will be a discrete set of modes
analogous to the biharmonic eigenfunctions discussed earlier in this
chapter, each of which will be stable at low γ̇ and attain a posi-
tive growth rate ω for sufficiently large shear rate. In the program
Figure_418_stretchcoil.m we find the different critical values Σ∗n
at which ω = 0 by appropriate choice of the initial guess for the
iterative solution.
n=1
-1
!:=2 0 :=2
@
Such analysis yields the critical value Σ∗ = −0.3932 for the first
unstable mode, Σ2∗ = −1.9876 for the second mode, Σ3∗ = −4.955 for
the third, etc. Those higher order shapes (Fig. 4.18) are visible in the
montage of filament shapes shown in Fig. 4.19. The amplitude of
the buckled filaments, taken as 1 − L/L, where L is the end-to-end
displacement, displays a thermal-fluctuation-rounded pitchfork bi- 20
H. Manikantan and D. Saintillan.
furcation (Fig. 4.19). A theoretical treatment of the combined effects Buckling transition of a semiflexible fil-
ament in extensional flow. Phys. Rev. E,
of thermal fluctuations and nonlinearities on the stretch-coil transi- 92:041002, 2015
tion has recently been developed.20
x At h ΓL2
χ= , τ= , H= , σ= , (4.106)
L ζ ⊥ L4 L A
The first two surface terms vanish identically at both ends by virtue
of the boundary conditions, as does the third term Hχ Hτ at χ = 0,
but Hχ (1, τ ) Hτ (1, τ ) is not forced to vanish by any boundary condi-
tions. Using (4.107) to simplify the integral in (4.109) we obtain
Z 1
Eτ = − dχ ( Hχχχχ + σHχχ)2 − σHχ (1, τ ) Hτ (1, τ ). (4.110)
0
While the integral is clearly nonpositive, the sign of the surface term
is indefinite, so we can not conclude that Eτ < 0 and indeed the
surface term represents the injection of energy into the filament that
sustains oscillations. Stated another way, the surface term breaks the
self-adjointness of the Sturm-Liouville operator ∂χχχχ + σ∂χχ , and
therefore its eigenvalues need not be real.
78 biological physics and fluid dynamics
0 Figure_421_and_423_followerforce.m.
Re(!)
-100
0 10 20 30 40 50
<
0 0.5 1
@
σ = 45, and plot it during one cycle of motion. As shown in Fig. 4.23
the dynamics have the character of flagellar motion.
-1
-2
0 0.5 1
@
5
Surfaces & Membranes (Lectures 13-16)
gαβ = xα · xβ , (5.1)
5.1.2 Catenoids
Now we ask how these geometric quantities arise in variational cal-
culations involving surfaces. The simplest context is that of surfaces
endowed with an energy per unit area (surface tension), for which
the area functional A[ζ ] is of interest,
Z d q
A[ζ ] = dz ζ 1 + ζ z2 . (5.12)
−d
δA
= −2ζ H, (5.13)
δζ
Viewing the radius of the hoops as a fixed quantity and the sepa-
ration 2d as adjustable, we see in Fig. 5.2(a) that this transcenden-
tal equation has two solutions for small D that merge into one at
D = Dc = 0.6627, and there are no solutions for larger D. The two
branches of solution, labelled α± , are shown in Fig. 5.2(b) to merge
together in what is known as a saddle-node bifurcation.
1
(a) (b)
,+
0.5
, cosh(D=,) ! 1
0.75
0:8
Dc ,c
0 0.5
,
0:6
0:5
(α− ), and is the true least-area surface. The areas of the two branches
of solution can be calculated analytically in terms of the solutions
α± . Expressed in dimensionless form as A± = A± /R2 , they are
α± 2D 2D
A± = sinh + , (5.18)
2 α± α±
84 biological physics and fluid dynamics
and are shown in Fig. 5.4. The upper and lower branches of solution
meet at the critical area Ac = 1.199 . . .. A more involved calculation
is necessary to show the stability of the two branches, which are
indicated in the figure. We shall return to these catenoidal solutions
in our study of membrane tethers in Sec. 5.3.2.
l e
0.75 ab
A
st
0.5
l D$ Dc
ca
r i ti
0.25 bc
su
0
0 0.2 0.4 0.6 0.8
D
R = R̄ + ∑ an sin qn z, (5.20)
n
1
R̄ ≃ R0 −
4R0 ∑ a2n . (5.22)
n
If we now compute the surface area to the same order and impose
volume conservation, we find that the surface energy γS is
π γL h i
γS = 2πLR0 γ +
2 R0 ∑ (qn R0 )2 − 1 a2n + · · · . (5.23)
n
The constraint of fixed volume has led to the crucial term −1 inside
the square brackets instead of simply a factor of q2 associated with
tension without conservation of volume. Its significance can be seen
if we apply equipartition to this energy, from which we conclude that
k B TR0
⟨ a2n ⟩ = h i. (5.24)
πLγ (qn R0 )2 − 1
then the constraint of fixed enclosed area yields the same result as in
(5.22), and the line energy γL is
πγ
2R0 ∑
2
γL = 2πR0 γ + n − 1 a2n , (5.26)
n
and hence
k B TR0
⟨ a2n ⟩ = . (5.27)
πγ (n2 − 1)
86 biological physics and fluid dynamics
1
E ≃ 8πk c + k c ∑ l (l − 1) (l + 1) (l + 2) | alm |2 . (5.30)
2 l,m
R
σ dS, and assuming small gradients, we have
2
1
Z
E≃ d2 x k c ∇2 h + σ (∇h)2 . (5.31)
2
2d
If we decompose the shape in the usual manner as
and by equipartition
kB T
⟨|ĥ(q)|2 ⟩ = . (5.34)
A (k c q4 + σq2 )
It follows that for the unconstrained membrane the variance of dis-
placements is
k T qdq
Z
⟨ h2 ⟩ = B (5.35)
2π k c q4 + σq2
√
If σ = 0 (a free membrane), using the cutoff qmin = π/ A, we find
kB T
⟨ h2 ⟩ = A. (5.36)
4π 3 k c
The fact that the variance grows with the area A of the membrane
led Helfrich and Servuss to argue that there will be a scale A∗ at
which the variance will be of order d2 , at which point the membrane
strongly interacts with the walls. The precise value of the prefactor
c in the statement ⟨h2 ⟩ = cd2 varies with the assumptions of the
calculation, so we leave it unspecified at the moment. Then,
kc
A∗ = 4π 3 c d2 . (5.37)
kB T
The next step involves the idea that each patch of membrane of
area A∗ acts as a quasi-independent "particle" confined in the gap
of size 2d, bouncing around in thermal equilibrium like a particle
in a confined one-dimensional ideal gas. The pressure exerted on
the wall can be deduced from the average force the wall experiences
from collisions of the constituent "particles" with some assume mass
m. That force is is just the ratio of the momentum transfer ∆p =
m∆v = 2mv to the average time interval ∆t between collisions with
the wall. Since the average velocity vanishes by symmetry, we use
mv ∼ m⟨v2 ⟩1/2 , and set ∆t = 4d/⟨v2 ⟩1/2 and thus F = k B T/2d.
Dividing by the area A∗ we obtain the pressure P = F/A∗ , or
( k B T )2 1
P∼ . (5.38)
8π 3 ck c d3
The energy per unit area is the integral of P with respect to d,
1 kB T kB T
Vrep (d) = . (5.39)
16π 3 c k c d2
88 biological physics and fluid dynamics
where ζ is the drag coefficient, k is the trap stiffness, and ξ (t) is the
random noise. We define the Fourier transform pairs
Z
x̃ (ω ) = dte−iωt x (t), (5.42a)
dω iωt
Z
x (t) = e x̃ (ω ). (5.42b)
2π
Recall that the choice of where the factor of 2π is located in these
relations is a matter of convention, or it can be split symmetrically
√
between the two as a factor of 1/ 2π. The only constraint is that the
transform of the inverse transform yields the function, or
Z ∞
dωeiωt = 2πδ(t), (5.43)
−∞
(iωζ + k) x̃ (ω ) = ξ̃ (ω ). (5.44)
2k B Tζ 2D
P(ω ) = = 2 . (5.45)
2 2
ω ζ +k 2 ω + τ −2
k B T −t/τ
⟨ x (t) x (0)⟩ = e . (5.48)
k
The celebrated Wiener-Khinchine theorem, derived first by Wiener10 10
Wiener, N. Generalized harmonic
and extended by Khinchine,11 states that the power spectrum and analysis. Acta Math., 55:117–258, 1930
11
Khintchine, A. Korrelationstheorie
the autocorrelation function are Fourier transform pairs. With the der stationären stochastischen Prozesse.
optical trap example, this is verified by direct calculation: Math. Ann., 109:604–615, 1934
Z ∞
dω 2D k T
eiωt = B e−t/τ . (5.49)
−∞ 2π ω 2 + τ −2 k
Third, we can define a "susceptibility" Υ of the system as the con-
stant of proportionality between the response of the system to an
90 biological physics and fluid dynamics
external force and the amplitude of that force. From (5.44) we have
Υ = (1/ζ )/(iω + τ −1 ) which leads to the result that the power spec-
trum is proportional to the imaginary part of the response function,
2k B T
P(ω ) = − ℑ {Υ} . (5.50)
ω
This is the generalised Fluctuation-Dissipation Theorem.
With these definitions, we may now summarize the findings of
Brochard and Lennon. Using a spectrum analyser they determined
from the microscopic observations the quantities shown in Fig. 5.10.
d3
J=− ∇ p, (5.57)
12µ ⊥
and the depth-averaged velocity ū = J/d, which obeys Darcy’s Law,
J d2
ū = =− ∇ p. (5.58)
d 12µ ⊥
Equation 5.56 captures the essence of the dynamics relevant to
the flicker phenomenon, which requires an equation of motion for
deformations of the membrane shape away from the flat equilibrium.
The relevant PDE for the shape evolution can be derived in three
steps. First, note that conservation of fluid volume implies
∂
(d + 2h) = −∇⊥ · J. (5.59)
∂t
Second, from the lubrication theory above we have (5.57) relating the
flux to gradients in the pressure. Finally, the pressure p at the inter-
face is related to the displacement field h( x ) via the usual functional
derivative. In the long-wavelength approximation, this yields
p = k c ∇4⊥ h. (5.60)
Assembling all the pieces, and dropping the subscript ⊥ for conve-
nience, the equation of motion for h is
k c d3 6
ht = ∇ h. (5.61)
24µ
This unusual form, with 6 derivatives, thus arises from a combina-
tion of a conservation law, lubrication flow, and elasticity. By dimen-
sional analysis we deduce immediately from (5.65) that for perturba-
tions at frequency ω there is a length scale
1/6
k c d3
ℓ(ω ) = , (5.62)
24µω
92 biological physics and fluid dynamics
h = h0 e−iq·x−iωt , (5.63)
indicating exponential decay of the modes for all q. They are stable.
To find the response function Υ, we introduce an external pressure
field pext that forces the system, and obtain the dynamics
k c d3 6 d3 2
ht − ∇ h= ∇ pext . (5.65)
24µ 24µ
d3 q2 1
Υ(q, ω ) = . (5.66)
24µ iω + kc d3 q6
24µ
k B Td3 q2 /12µ
⟨|δd(q, ω )|2 ⟩ = 3 2 . (5.67)
ω 2 + k24µ
cd
q12
d2 q
Z
C ( R, ω ) = ⟨|δd(q, ω )2 ⟩eiq·R . (5.68)
(2π )2
and the prefactor sets G (0) = 1. Figure 5.11 shows the scaling func-
tion G (ξ ) obtained by numerical integration of (5.74). This result con-
Figure_511_flicker.m.
0 2: 4:
@ = R=`(!)
firms that the correlation function depends only on the ratio R/ℓ(ω ),
as found experimentally, and that therefore ℓ(ω ) is the length scale
governing the decay of correlations. The prefactor ℓ2 (ω )/ω displays
the ω −4/3 dependence seen in experiment, and the scaling function
G (χ) shows the exact experimental form shown in Fig. 5.10(b).
γR20 ΠR30 s
K = κR0 , T= , P= , η= , (5.76)
A A R0
1
Kηη + K3 − TK − P = 0, (5.77)
2
or, equivalently,
1
θηηη + θη3 − Tθη − P = 0. (5.78)
2
Given the curvature as a solution to (5.77), we obtain the θ via
Z η
θ (η ) = K ( η ′ ), (5.79)
0
1
θ (0) (η ) = η, T (0) = − P (0) , (5.82)
2
Lθ (1) (η ) = 0, (5.83)
surfaces & membranes (lectures 13-16) 95
3
Lθ (2) = P(2) + T (2) − m2 cos2 (ms), (5.86)
2
(1) 2
where the final term arises from a contribution (3/2)θη . Here we
must be careful to recognise that since cos2 (ms) = (1 + cos(2ms))/2,
the r.h.s. of (5.86) has a contribution in the nullspace of the operator
L, which involves the functions θ (η ) = 1, cos(ms), sin(ms). For there
to be a bounded solution at this order, the coefficients of all such
terms on the r.h.s. must vanish. Here, that condition requires the
constant term to vanish, and thus
3 2
P (2) + T (2) = m . (5.87)
4
1
θ (2) ( η ) = sin(2mη ). (5.88)
8m
Using the previously obtained results and now requiring the absence
of a term proportional to cos(ms) we find
3 2
T (2) = m +1 , (5.90)
8
ϵ2
θ (η ) = ϵ sin(ms) + sin(2ms) (5.91)
8m
3 2:
1
2
K(2)
3(2)
1 : 0
Y
0
-1
-1 0 -1 0 1
0 : 2: 0 : 2: X
2 2
Figure 5.12: Shapes of two-dimensional
vesicles under pressure. Panels show
(a) scaled curvature K (η ), (b) tangent
5.3.2 Tether formation angle θ (η ), and (c) shape of vesicle
for three values of ϵ (ϵ = 0 (black),
A classic problem in the biological physics of membranes is the for- 0.5 (blue) and 1 (red). Obtained with
Figure_512_2dvesicles.m.
mation of "tethers" under the action of a localized force. There are
many examples of this phenomenon in living cells, typically associ-
14
Fygenson, D.K., Marko, J.F. and
ated with cytoskeletal filaments pushing against the cell membrane. Libchaber, A. Mechanics of micro-
A controlled example of this was achieved some time ago14 using an tubule-based membrane extension.
optical trap to pull a tether from a vesicle (Fig. 5.13). In what fol- Phys. Rev. Lett., 79:4497–4500, 1997
z−c
ζ (z) = a cosh . (5.92)
a
D
0:6
for each chosen loop ratio β the parameter α necessary to match the 0:2
0:4
boundary conditions varies in the range (0, β) along the lower branch 0: 0
5
1
Z Z
E= kc dS(2H )2 + σ dS (5.95)
2
for an axisymmetric surface with boundary conditions ζ (0) = R,
ζ (d) = r0 , H (0) = 0 and the tangent plane at z = d is parallel to
the plane. The Euler-Lagrange equation for the Helfrich functional
is an amusing exercise in differential geometry which we leave to
the student. The result can be expressed in a compact form that has
obvious similarities to the equivalent result for planar vesicles (5.75),
k c ∇2 H + 2H 3 − 2HK − σH + p = 0, (5.96)
where p is the pressure difference between the the inside and outside
of the tether. Here ∇2 is the covariant Laplacian. We shall ignore any
pressure difference in this problem, as it can be shown to introduce
only minor corrections to the main results.
98 biological physics and fluid dynamics
where
kc
ϵ= . (5.98)
σR2
We shall be interested in the large-tension limit ϵ ≪ 1, as this gives
narrow tethers. Written this way, with the small parameter multiply-
ing the highest derivative in the problem, it is clear that a boundary
layer is likely to form, as the introduction of bending energy is a
singular perturbation. The limit ϵ = 0 is just the minimal surface
problem, so the choice H = 0 at the bottom loop is consistent with
that limit, and no boundary layer will form there. But the tangent
condition at the upper loop is inconsistent with H = 0 and a bound-
ary layer will form there. Such a layer will be such that the elastic
and surface tension contributions balance (ϵ∇2 H ∼ H, giving a (di-
√ √
mensionless) length scale ∼ ϵ (dimensionally, R ϵ).
To understand the boundary layer problem, we first work in the
Monge representation that is suitable for shallow surfaces slightly
deformed from the disc at z = 0, so D ≪ 1. A useful pair of coor-
dinates (ξ 1 , ξ 2 ) to describe the surface involves the arclength s along
the meridian and the angle φ around the z-axis, so that points on the
surface are given by
With the previous convention for raising and lowering indices and
gik gkj = δji one finds the mean curvature
1 ij 1 rss zs
H= g Kij = − , (5.102)
2 2 zs r
Gaussian curvature
rss
K=− , (5.103)
r
and Laplacian
1 √ 1 d d
∇2 = √ ∂i gij g∂ j = r . (5.104)
g r ds ds
Thus, if the surface is described by the height function z(r ) above the
plane z = 0, then to leading order s ∼ r, H ∼ (1/2)∇2 z, and the
shape equation becomes
ϵ∇4 z − ∇2 z = 0, (5.105)
where
21 d d
∇ = r . (5.106)
r dr dr
surfaces & membranes (lectures 13-16) 99
zout (r ) = b1 + b2 ln r. (5.107)
This is solved by
zin
= c1 + c2 ln ρ + c3 I0 (ρ) + c4 K0 (ρ), (5.110)
d
where I0 and K0 are modified Bessel functions. The coefficient c3
vanishes due to the divergence of I0 as ρ → 0. The boundary condi-
tions at r = r0 yield two constraints, from which we deduce
zin ρ
= 1 + c4 K0 (ρ) − K0 (ρ0 ) + ρ0 K1 (ρ0 ) ln , (5.111)
d ρ0
√
where ρ0 = r0 / ϵ. Matching zin and zout , the constant terms of zin
must vanish and the coefficient of the log term matches b2 , yielding
zout ρ0 K1 (ρ0 ) ln r
= (5.112)
d ρ0 K1 (ρ0 ) ln r0 + K0 (ρ0 )
and
zin ρ0 K1 (ρ0 ) ln r + K0 (ρ)
= (5.113)
d ρ0 K1 (ρ0 ) ln r0 + K0 (ρ0 ) .
0
-1 -0.5 0 0.5 1
r
100 biological physics and fluid dynamics
0.5
z
-1 -0.5 0 0.5 1
r(z)
Cell motion can take many forms, from crawling on solid substrates
by means of membrane extensions ("filopodia" or "lamellipodia"), mi-
gration within confluent tissues, to swimming through fluids through
the action of moving appendages. Within this course we focus on
swimming, as it provides a context that is rich in phenomenology
and has proven to be amenable to substantial quantitative analysis.
250 μm
frequencies that can reach 102 − 103 s−1 by the action of rotary mo-
tors embedded in the cell wall, powered by a flux of proteins passing
through them. Eukaryotes, as exemplified by uni- and multicellular
algae, choanoflagellates, etc. achieve motility through the undula-
tions of flagella driven by the action of motor proteins that slide the
underlying filaments within the flagellum past each other through
conformational changes powered by ATP consumption. These wave-
forms consist typically of a "power stroke", with a nearly straight
flagellum sweeping through a large angle, followed by a "recovery
stroke", in which a tightly curled flagellum is brought back to its
starting point by moving close to the cell body. Both large unicellular
organisms such as Paramecium and certain multicellular organisms
composed of tens, hundreds, or even thousands of flagellated cells
are often termed "ciliates", and the collective beating of their flag-
ella may exhibit perfect frequency and phase locking, or frequency
locking and phase modulations termed "metachronal waves".
The regime of Reynolds numbers inhabited by various swimming
organisms is summarized in Table 6.1; there is a clear distinction
between the Stokesian regime of the flagellated organisms in Fig. 6.1
and that of the larger swimmers familiar to us on a human scale.
In the absence of any body forces and torques acting on the swim-
mer, the unknowns are found by demanding that the total force and
torque on the swimmer vanish,
ZZ
F( t ) = dS σ · n̂ = 0, (6.3a)
Z ZS
T( t ) = dS r × (σ · n̂) = 0, (6.3b)
S
cellular motion (lectures 17-20) 105
for some tensor H(r) that does not depend explicitly on time and is
found as a solution to the swimming problem.
The total distance travelled is
Z t
1
∆= dt U(t), (6.8)
t0
This is illustrated in Fig. 6.4. Neither (a) a pair of spheres that cycli-
cally move together and apart, nor (b) a hinged object that cyclically
(c)
opens and closes (i.e. a scallop) can achieve net locomotion. But, (c)
a three-link swimmer can swim by following a cycle of shapes.
The need to break time-reversal invariance to achieve locomotion
explains the ubiquitous appearance of helical prokaryotic flagella ro-
tated along their long axis or undulating eukaryotic flagella with
either sperm-like sinusoidal traveling waveforms (as we saw in our Figure 6.4: Swimming sequences and
the scallop theorem.
study of wiggling elastica in Sec. 4.3.2), or two distinct phases of
beating with very different waveforms. Two classic examples of the
latter is shown in Fig. 6.5, taken from the very early study by James
Gray of the ciliated gills of the blue mussel Mytilus edulis6 and very 6
J. Gray. The mechanism of ciliary
movement. Proc. R. Soc. Lond. B, 93:104–
recent work on the green alga Chlamydomonas.7 The first phase of
121, 1922
a beat cycle is the "power stroke", in which a nearly straight cilium 7
Leptos, K.C. and Chioccioli, M. and
roughly pivots around its attachment point to the underlying tis- Furlan, S. and Pesci, A.I. and Gold-
stein, R.E. Phototaxis of Chlamydomonas
sue or cell, followed by a "recovery stroke" in which the waveform is
arises from a tuned adaptive pho-
highly curled and returns to the original position close to the surface. toresponse shared with multicellular
Volvocine green algae. Phys. Rev. E, 107:
014404, 2023
Figure 6.5: Ciliary waveforms. Left
panel: sketch of the waveform in the
mussel Mytilus edulis. (a) Power stroke.
(b) Recovery stroke. From Ref. 6.
Right panel: flagellar waveforms of
the two flagella of Chlamydomonas rein-
hardtii. Cis (red) and trans (blue) flagella
are shown during one complete cycle.
From Ref. 7.
∇4 ψ = 0. (6.13)
∂ψ( xs , ys )
u x ( xs , ys ) = = 0, (6.14a)
∂y
∂ψ( xs , ys ) ∂ys
uy ( xs , ys ) = − = = −ϵ cos ϕ, (6.14b)
∂x ∂t
∂ψ( x, y → ∞)
u x ( x, y → ∞) = = U, (6.14c)
∂y
∂ψ(y → ∞)
uy ( x, y → ∞) = − = 0, (6.14d)
∂x
108 biological physics and fluid dynamics
and note that ∇4 ψ(n) = 0 ∀n. When expanding the boundary con-
ditions evaluated on the sheet we make use of relationships familiar
from the boundary perturbation theory used in Sec. 2.4.2, such as
∂ψ(1)
ψ(1) ( xs , ys ) = ψ(1) ( xs , 0) + ϵ sin ϕ ( x s , 0) + · · · . (6.16)
∂y
∂ψ(1) ∂ψ(1)
( x, 0) = 0, ( x, 0) = cos ϕ,
∂y ∂x
∂ψ(1) ∂ψ(1)
( x, ∞) = 0, ( x, 0) = U (1) . (6.17)
∂x ∂y
∂ψ(2) ∂ 2 ψ (1)
( x, 0) = − sin ϕ ( x, 0), (6.20)
∂y ∂y2
which are proportional to sin2 ϕ. This and the other boundary con-
dition at ( x, 0) force the system at its first harmonic, leading to
(and denote the average by angular brackets), the periodic terms give
no contribution and thus
(2)
(2) ∂ψ
U = , (6.22)
∂y
∂ψ(2) ∂ 2 ψ (1)
U (2) = ( x, 0) = − sin ϕ ( x, 0 ) . (6.23)
∂y ∂y2
tion to the case of two nearby sheets separated by a mean distance Figure 6.7: Two nearby waving sheets
separated on average by 2h and with a
2h and offset by a phase angle 2Θ, as in Fig. 6.7. Here we quote the
phase shift Θ.
essential results, leaving algebraic details to the reader.
Having rescaled x by the wavevector k in (6.10), if the two sheets
are separated on average by 2h, then we define the dimensionless
δ = kh, (6.25)
y = δ + y1 = δ + ϵ sin (ϕ + Θ) ,
y = −δ + y2 = −δ + ϵ sin (ϕ − Θ) , (6.26)
(1) ∂y1
Ẇ (2) = − p1 . (6.34)
∂t
If we user an overbar to denote an average over a temporal cycle
of oscillations, and define the (scaled) mean rate of working W as
µkW (2) = Ẇ (2) , then we obtain
Since it is evident from (6.33) that β > α, ∀δ, we deduce that W (Θ, δ)
is minimised for Θ = 0; the in-phase synchronised state has the
lowest dissipation rate.
Intuitively, since the dissipation rate is a spatial integral of the
squares of fluid velocity gradients, it stands to reason that the con-
figuration in which the separation between the sheets is nearly con-
stant throughout the waveform will be less dissipative than states 1
with large thickness variations. This can be made concrete by con-
sidering the ratio R(δ) = W (0, δ)/W (π/2, δ)) of dissipation rates
for in-phase and out-of-phase configurations. We find
R
0.5
sinh δ cosh δ − δ
R(δ) = tanh2 δ , (6.36)
sinh δ cosh δ + δ
which is plotted in Fig. 6.8. This ratio tends to unity at large separa- 0
0 1 2 3 4
tion, indicating in that limit that the two sheets do not interact and / = kh
the dissipation is independent of the phase shift. But for δ < 2 (i.e. Figure 6.8: Dissipation ratio in Tay-
lor’s calculation, as a function of
h/λ < 1/π ∼ 0.3), R drops to small values, signalling a significant scaled sheet separation. Obtained with
streamlining effect associated with synchronisation. Figure_608_Taylordissipation.m.
cellular motion (lectures 17-20) 111
motor proteins embedded in the cell wall. For eukaryotes like sper-
matozoa and choanoflagellates, the flagellum has a sinusoidal shape
deformed by the sliding action of motor proteins along the length of
the filament. As shown in Fig. 6.10, in both cases we find segments
of the flagellum that are moving down when tilted upward from left
to right, and moving up when tilted downward from left to right.
Consider the case shown in Fig. 6.11, and project the upward
speed U onto the tangent and normal directions, giving U∥ and U⊥ .
The drag forces F∥ and F⊥ are as indicated. The vectorial sum of par-
Figure 6.11: Normal and tangential
allel and perpendicular forces necessarily has a component (shown forces on a moving filament segment.
in green) that is to the left, opposite to the direction of the propaga-
tion of the traveling wave of deformation along the filament.
We can verify that this conclusion holds for any inclination angle
of the segment by considering a straight segment oriented as in Fig.
6.12. Here, the tangent vector is t̂ = cos θ, sin θ êy and the segment
velocity is u = (0, u). A simple calculation yields the force
f = ζ ⊥ − ζ ∥ u sin θ cos θ êx
h i
+ −ζ ⊥ u + ζ ⊥ − ζ ∥ u sin2 θ êy , (6.41) Figure 6.12: A tilted rod segment mov-
ing upwards.
and thus f · êx = ζ ⊥ − ζ ∥ u sin θ cos θ. We deduce that if u > 0
and sin θ < 0 we obtain r · êx < 0 and if u < 0, sin θ > 0 and again
r · êx < 0; the swimming force component is independent of the sign
of u and propulsion is directed forward everywhere. For a helix
every point is the same... for a waving flagellum they are different,
but all contribute in the same sense.
are inclined at the same angle θ with respect to the x −axis, we can
use the results derived in the context of Fig. 6.12, where the upward
speed of the segment shown is u = −ΩR, to deduce that the total
drag force acting along x is
F = f x L = − ζ ⊥ − ζ ∥ sin θ cos θRLΩ, (6.47)
where
A = ζ ⊥ sin2 θ + ζ ∥ cos2 θ L, (6.51a)
B = ζ ⊥ − ζ ∥ sin θ cos θRL, (6.51b)
D = ζ ⊥ cos2 θ + ζ ∥ sin2 θ R2 L. (6.51c)
cellular motion (lectures 17-20) 115
or, equivalently,
! ! !
A0 + A B U Bω
=− (6.54)
B D0 + D Ω Dω.
Figure 6.16 displays these results for the case θ = π/6 and various
0.2 1 0
~ !1
(a) (b) (c)
body rotation frequency +
~
;c = 2
;c = 0:5
0 0 -1
0 10 20 0 10 20 0 10 20
` ` `
Figure 6.16: Swimming dynamics of
a model bacterium, from Eqs. (6.57).
values of ρc . For each cell body size there is an optimum flagel- Panels display the rescaled (a) swim-
lar length for maximum swimming speed, followed by a slow decay ming speed, (b) cell body rotation fre-
quency, and (c) flagellar rotation fre-
in the speed with increasing length. As ℓ increases, the body ro- quency in lab frame. Obtained with
tation frequency tends toward the negative of the flagellar rotation Figure_616_bacterialswimming.m.
2a
Figure 6.17: Two coupled bead-spring
F oscillators in the NEL model. Each
bead (red circles) is pushed around an
orbit by a constant tangential force F
R1 and is held by a spring (green) to the
φ� φ�
R2 equilibrium orbit radius R0 .
F
R0 R0
x l
and likewise for particle 2. Here, we take the flow to be the far-field
flow due to a moving sphere,
s · (I + n̂12 n̂12 )
v12 = +··· (6.65)
|r12 |
where the vector r12 points from 2 → 1 and
3a 3a
s= ṙ ≃ R2 φ̇2 ê φ2 , (6.66)
4 12 4
the latter approximation holding under the assumption of small dis-
placements from a circular orbit.
In the "weak-coupling" regime ℓ ≫ R the terms simplify:
where
R0 R2
φ̇1 = ω1 − ρJ ( φ1 , φ2 ) φ̇2 , (6.73)
R1 R1
Ṙ1 = −τλ−1 ( R1 − R0 ) + ρR2 K ( φ1 , φ2 ) φ̇2 , (6.74)
This and the equivalent relation for R2 can now be substituted into
the angular dynamics and the result simplified under the assumption
cellular motion (lectures 17-20) 119
from which we infer that the synchronized state is an attractor for all
initial conditions. From the form of γ we see that the time to achieve
synchrony can be significantly longer than an orbital period if the
oscillators are far apart.
lags leads
speeds up
cell body (Figs. 6.21(a-l)). In the simplest analysis, the peaks in inten-
sity as flagella move through the regions define the endpoint of each
cycle, with the phase angle linearly interpolated in time between.
The data in Fig. 6.21(m) on the scaled phase difference ∆(t) =
( φ1 (t) − φ2 (t))/2π shown in blue reveals that within periods of syn-
chrony there are significant fluctuations away from the mean value
of ∆ (red). A simple calculation shows that the fluctuations observed
are too large to arise from timing jitter due to purely equilibrium
thermal fluctuations of elastic flagella. The data also show a transi-
tion over the course of ∼ 1 s of the running mean value of ∆ by one
unit, followed by resynchronisation at that new value. This transition
corresponds to the loss of one complete cycle of motion for flagellum
2 relative to flagellum 1. In the language of coupled oscillators this is
termed a "phase slip", and corresponds to a transition from one min-
imum to an adjacent one of the effective potential V (θ ). This kind of
cellular motion (lectures 17-20) 121
dV (θ )
θT = − + ξ ( t ), (6.83)
dθ
where the noise has the typical properties we discussed in (3.45) of
Sec. 3.3, ⟨ξ (t)⟩ = 0 and ⟨ξ (t)ξ (t′ )⟩ = 2k B Teff δ(t − t′ ), where Teff is
some effective temperature. This stochastic Adler equation allows
for a heuristic view of synchronisation, in that we can view (6.83)
as the overdamped equation of motion of a particle with coordinate
θ on the "landscape" defined by V (θ ), subject to noise. Periods of
oscillator phase synchrony then correspond to sojourns in one of
the minima of V, and the noise induces quasi-thermal fluctuations
in a harmonic potential (like the optical trap considered in (5.41)
of Sec. 5.2.3 on flicker phenomena of membranes). The linearised
dynamics (6.80) around the stable minimum θ− ∗ implies that there is
k B Teff ′ 1
⟨ψ(t)ψ(t′ )⟩ = p e−|t−t |/τ , τ= p . (6.85)
γ2 − (δω )2 γ2 − (δω )2
Phase slips like that shown in the middle panel of Fig. 6.21 are
thermally-assisted hops from one minimum of V to an adjacent one,
122 biological physics and fluid dynamics
R± = e2πδ/Teff . (6.86)
Thus, measurement of (i) the amplitude, (ii) the decay time of the
autocorrelation function, and (iii) the relative probability of left-right
hopping, provides three relations involving the three parameters of
the model (γ, δ, Teff ), allowing all to be determined from the data.
Finally, we note that the NEL model makes testable predictions
regarding the coupling strength due to hydrodynamic interactions.
For example, it shows that the coupling constant should vary with 18
Goldstein, R.E., Polin, M. and Tuval,
flagellar length as γ ∼ L3 ; experiments that take advantage of the I. Emergence of Synchronized Beat-
ing during the Regrowth of Eukaryotic
ability of Chlamydomonas to shed and then regrow flagella have been Flagella. Phys. Rev. Lett., 107:148103,
used to confirm this prediction.18 Second, the dependence of the 2011
coupling on the oscillator separation, γ ∼ 1/ℓ, a consequence of the 19
Brumley, D.R., Wan, K.Y., Polin, M.
hydrodynamical coupling due to stokeslet flows, has been confirmed and Goldstein, R.E. Flagellar synchro-
nization through direct hydrodynamic
by examining the interaction between flagella on two separated cells interactions. eLife, 3:e02750, 2014
interacting only through the fluid between them.19
6.4.1 Stokeslets
We begin by recalling the flow field induced by a point force, u(x) =
G(x − x′ ) · F, with the Oseen tensor
1
G(r) = (I + r̂r̂) , (6.87)
8πµr
where r = |r| and r̂ = r/r is a unit vector. We refer to this flow field
as a "stokeslet". In terms of components, we have
Fj ri r j
δij
ui = + 3 , (6.88)
8πµ r r
so for a force along z, the flow in the x − z plane has components
Figure 6.22: Streamlines of a stokeslet
u x = ( F/8πu) xz/r3 and uz = ( F/8πµ)(1/r + z2 /r3 ), giving rise to
at the origin, pointing upwards in
the streamlines shown in Fig. 6.22. the x − z plane. Obtained with
A biophysical realisation of this flow is seen with an organism like Figure_622_stokeslet.m.
Volvox in Figs. 1.2 and 6.1, as its density ρV is slightly greater than
ρw , that of the water it swims through. It is therefore acted on by a
downward (negatively) buoyant force
4π
Fb = − δρw gR3 ez , (6.89)
3
cellular motion (lectures 17-20) 123
antisym P (d × e) × r
udipole (r; e, d) = , (6.94a)
8πµ r3
P 3(e · r)(d · r)r (d · e)r
sym
udipole (r; e, d) = + . (6.94b)
8πµ r5 r3
a3
u= Ω × r, and T = 8πµa3 Ω. (6.95)
r3
If we substitute the second relation into the first to eliminate refer-
ence to the sphere, we obtain the "rotlet" flow uR (r; T), defined as the
flow due to a point torque T,
1 T×r
uR (r; T) = , with T = P d × e. (6.96)
8πµ r3
cellular motion (lectures 17-20) 125
6.26 shows the flows fields for pushers and pullers. The streamline
geometry is the same in the two cases, only the direction of motion
along them switches. A characteristic feature of the flows is the nodal
lines at angles θn such that cos2 θn = 1/3 or θn ≃ ±54.7◦ .
These results can be compared to experimental measurements of 22
Drescher, K., Dunkel, J., Cisneros,
L.H., Ganguly, S. and Goldstein, R.E.
the flow field around individual freely-swimming E. coli bacteria Fluid dynamics and noise in bacte-
shown in Fig. 6.27.22 These data were obtained by tracking the rial cell-cell and cell-surface scattering.
motion of micron sized tracer spheres around individual cells swim- Proc. Natl. Acad. Sci. USA, 108:10940–
10945, 2011
ming in the focal plane of the microscope to obtain short traces of
the fluid velocity field, and then averaging over a very large number
of traces rotated to align the cellular axis in a common direction. The
126 biological physics and fluid dynamics
ṙ = Up + u(r), (6.103)
128 biological physics and fluid dynamics
where u(r ) is the flow evaluated at the centre of mass of the particle.
The subtlety in this approach is that the presence of the particle will
disturb the flow, so even in the simple situation where we specify
a flow field u∞ (r) far away from the cell it is unclear the level of
approximation entailed in setting u(r) = u∞ (r).
A systematic approach to this is given by Faxén’s laws,25 whose 25
Faxén, H. Der Widerstand gegen
derivation we summarise below (see also Ref. 1, Sec. 10.1). We die Bewegung einer starren Kugel
in einer zähen Flüssigkeit, die zwis-
assume that the flow u∞ is that specified in the absence of the swim- chen zwei parallelen ebenen Wänden
mer, and that in the absence of the flow a rigid spherical swimmer eingeschlossen ist. Ann. Phys., 373:89–
119, 1922
moves with velocity Us and angular velocity Ωs For a rigid swimmer,
the principle of superposition is valid, so that if U f and Ω f are the
velocities due to the flow, then the total velocities are
U = Us + U f , (6.104a)
Ω = Ωs + Ω f . (6.104b)
3µ
σ̂ · n = − Û − 3µΩ̂ × n, (6.108)
2a
1
Z
Uf = dSu∞ , (6.109a)
4πa2
3
Z
Ωf = dSn × u∞ . (6.109b)
8πa3
These very powerful results express the kinematics of a rigid moving
sphere in terms of the imposed flow (not the actual flow), evaluated
on its surface. Importantly, that imposed flow is defined everywhere
in space, including inside the sphere.
Recalling that we wish to relate the motion of the sphere to the
value of the fluid velocity at its centre, while (6.109a) requires it on
the sphere surface, we Taylor expand the velocity around the centre.
By symmetry the relevant terms involve higher and higher order
contributions of the form ∇2n u∞ , which all vanish for n ≥ 2 by the
Stokes equation, leaving a simple relation between the translational
velocity and the external flow evaluated at the particle position. This
leads leads to the positional dynamics
a2 2
ṙ = Us + u∞ (r) + ∇ u∞ (r ). (6.110)
6
1
Ωf = ∇ × u∞ (r ). (6.111)
2
We conclude that the evolution equation for the orientation vector of
a spherical particle, the partner relation to (6.110), is
1
ṗ = (∇ × u∞ (r)) × p. (6.112)
2
Already in our discussion of the instability of semiflexible fila-
ments under compressional flows (Fig. 4.19) we saw that elongated
objects can rotate in flows that lack vorticity but instead have a finite
rate-of-strain. This implies that there must be an additional contri-
bution to the evolution equation for p beyond that in (6.112). It also
raises the question of whether for elongated self-propelled swim-
26
Hohenegger, C. and Shelley, M. J. Dy-
mers there is an additional contribution beyond (6.110). A careful namics of complex biofluids. In New
analysis26 of these questions is summarized below. Trends in the Physics and Mechanics of Bi-
To begin, we define xc (t) to be the position of the centre of a rod ological Systems: Lecture Notes of the Les
Houches Summer School: Volume 92, July
of length L, parameterized by arclength s ∈ (− L/2, L/2) relative to 2009. Oxford University Press, 2011
that centre. Choosing the orientation vector p to lie along the rod’s
130 biological physics and fluid dynamics
u∞ (r) = A · r, (6.115)
U f = ẋc = A · xc , (6.118)
L3
Z L/2
ds sp(t) × sf̃(t) = p(t) × f̃(t) = 0, (6.119)
− L/2 12
which shows that f̃ ∥ p, and we may therefore write f̃(t) = Γ(t)p(t)
for some function Γ.
The calculation is completed by noting that the relation (6.116) can
be inverted (as we did in deriving (4.99)) to obtain
1
urod − u∞ = − (I + pp) · f, (6.120)
ζ⊥
and thus, with f simplified as above, we must have
sΓ 2sΓ
s (ṗ − A · p) = − (I + pp) · p = − p. (6.121)
ζ⊥ ζ⊥
Hence, ṗ − A · p = −(2Γ/ζ ⊥ )p. Taking the dot product with p and
noting that p · ṗ = 0 since p is a unit vector, we obtain
ζ⊥
Γ= p · A · p, (6.122)
2
cellular motion (lectures 17-20) 131
and finally
ṗ = (I − pp) · (A · p) . (6.123)
With some further vector calculus manipulations one can show that
A can be expressed in terms of the rate-of-strain tensor E and vor-
ticity ω such that the final evolution equation for p takes the form
1
ṗ = Ω f × p, Ωf = ω + p × (E · p) . (6.124)
2
The result (6.124) has the very same vorticity contribution as the rigid
sphere and an additional contribution from the rate-of-strain tensor.
This second contribution describes the alignment effect we saw in
experiments on actin filaments discussed in Chapter 4.
Thus far, we have obtained the contribution to the angular velocity
from the rate-of-strain tensor for spheres, where there is no contri-
bution, and for extremely elongated rods, where we found (6.124).
These results are special cases of a more general result due to Jef-
fery,27 who computed the contribution for spheroids Combining his 27
Jeffery, G.B. The motion of ellipsoidal
result with bottom-heaviness and the vorticity contribution, we ob- particles immersed in a viscous fluid.
Proc. R. Soc. Lond. A, 102:161–179, 1922
tain "Jeffery’s equation",
1 1
ṗ = p × (ez × p) + ω × p + αp · E · (I − pp) , (6.125)
τbh 2
where, α is a parameter of the ellipsoid given in terms of the ratio r
of the major to minor axes,
r2 − 1
α= , (6.126)
r2 + 1
with α = 0 for a sphere and α → 1 for a highly elongated object.
There are many interesting and biologically relevant results that
follow from Jeffery’s equation (6.125) applied to a range of flow ge-
ometries and particle shapes. In the following two subsections we
illustrate a few key results for swimmers far from surfaces.
ẋ = U cos ϕ, (6.129a)
x2
ż = U sin ϕ + U0 1 − , (6.129b)
L2
xU0 1
ϕ̇ = − + cos ϕ. (6.129c)
L2 τbh
132 biological physics and fluid dynamics
x z t
X= , Z= , T= , (6.130)
L L τbh
leading to
XT = β cos ϕ, (6.131a)
ZT = β sin ϕ + β 0 1 − X 2 , (6.131b)
ϕT = − β 0 X + cos ϕ, (6.131c)
where
U0 τbh Uτbh
β0 = , and β= (6.132)
L L
are the ratios of the bottom-heaviness relaxation time to the time
necessary for the maximum flow or the swimming organism to move
the channel half-width. Note that β > 0, while β 0 can be of either
sign; β 0 > 0 is an upward channel flow, β 0 < 0 is downward.
The dynamical system (6.131) allows us to understand the phe-
nomenon of gyrotaxis, the tendency of bottom-heavy organisms to
swim into regions of vorticity. If β 0 < 0 (downward flow), then we
can see from (6.131c) that ϕ will tend to increase when X > 0, mean-
ing that the cell rotates to point toward the flow centreline. Like-
wise, ϕ decreases when X < 0 leading again to motion toward the
center of the parabolic flow. For a flow with constant vorticity it is
possible for there to be a fixed inclination angle ϕ, but in the case
of Poiseuille flow the trajectories eventually reach the centreline with
the cell pointing upwards, but being advected downward at the same
time (if β 0 + β < 0), as depicted in Fig. 6.30. This leads to the phe-
nomenon of gyrotactic focusing of algae, first described by Kessler,
in which a thin concentrated plume of algae forms at the centreline.
:
(a) (b) Figure 6.30: Gyrotaxis of a bottom-
heavy organism in channel flow. (a)
2 Trajectory for downward parabolic
flow, showing motion to centreline,
with β = 2.5, β 0 = −3, starting at
X = 0.7, Z = 0, and ϕ = π/4 (red) and
X = −0.5, Z = 0, and ϕ = 5π/8 (blue).
(b) Corresponding tilt angles. Obtained
0 :=2 with Figure_630_gyrotaxis.m.
Z
-2
0
-1 0 -1 0 5 10 15
X T
cellular motion (lectures 17-20) 133
f˙ γ̇
=− , (6.137)
1+ f 2 r + r −1
0 1
Figure 6.31: Dynamics of the orienta-
tion angle in Jeffery orbits for various
aspect ratios of the cell. Curves are an-
gle ϕ in (6.138) for r = 1 (black), r = 2
(blue) and r = 10 (red). Obtained with
cos ?
-2: 0 Figure_631_and_632_Jefferyorbits.m.
?
(a) (b)
-4: -1
0 1 2 0 1 2
t=Torbit t=Torbit
rapid transitions to the next plateau. This is seen even more clearly
in Fig. 6.31(b) where we plot the projection cos ϕ. Note the important
feature that ϕ(t) does not depend on the vertical position y(t), sd the
vorticity and rate-of-strain are constant for this linear shear profile.
The swimming trajectories in this shear flow follow from the gen-
eral equation of motion (6.103) with the given ϕ(t). If we introduce
the scaled time T = t/Torbit and scale lengths with the quantity
ℓ = UTorbit , giving X = x/ℓ and Y = y/ℓ, then
XT = 2π (r + r −1 )Y + cos Φ( T ), (6.140a)
YT = sin Φ( T ), (6.140b)
0.4
Figure 6.32: Jeffery orbits in the
laboratory frame for r = 5. Curves
0.2
correspond to different initial po-
sitions as denoted by the solid
circles. Vertical scale is exagger-
Y
-0.2
-5 0 5 10 15
X
h
6.6 Surface-Mediated Interactions
Suppose first that the force points upward, away from the surface, so
Fj = Fδj3 . If the lateral position of the singularity is (xs = ys = 0),
and it is at zs = h, the (i = 1, 3) components of the flow are
x (z − h)
F x (z + h)
u⊥x = − , (6.144a)
8πµ r3 R3
( z − h )2 ( z + h )2
⊥ F 1 1
uz = − + − . (6.144b)
8πµ r R r3 R3
It is easily seen that uz vanishes identically at z = 0, where r = R,
as does ∂u x /∂z. Figure 6.34(a) shows the fluid velocity vector field
and heat map of velocity magnitude. While the fluid velocity falls
( z − h )2 ( z + h )2
∥ F 1 1
ux = + + + , (6.145a)
8πµ r R r3 R3
x (z − h)
∥ F x (z + h)
uz = + . (6.145b)
8πµ r3 R3
the streamlines of which are shown in Fig. 6.34(c). In this case, the
∥
parallel flow u x far above the singularity retains its 1/z falloff.
136 biological physics and fluid dynamics
where α = 1, 2. The quantity δjα δαk − δj3 δ3k is nonzero only when
j = k and is +1 for j = 1 or 2 and −1 for j = 3. The solution
involves terms with prefactors proportional to various powers of h.
That with no prefactor is just the image stokeslet as in (6.142). The
term proportional to 2hF has the interpretation of a force dipole (as it
is the derivative of terms associated with a force), while that scaling
as 2h2 F is a source dipole (the derivative of a source).
For the same vertically-oriented force (Fj = Fδj3 ), the flow is
x (z − h) x (z − h) 6hxz(z + h)
F
u⊥
x = 3
− − , (6.147a)
8πµ r R3 R5
(z − h)2 z2 + h2 6hz(z + h)2
F 1 1
u⊥
z = − + − − . (6.147b)
8πµ r R r3 R3 R5
This flow field is shown in Fig. 6.34(c), where we observe the new
feature of closed streamlines. These do not appear when the force is
parallel to the surface. For the perpendicular case, the existence of
recirculating flows is thought to have implications for filter-feeding
of organism that attach themselves to surfaces with a slender stalk.30 30
Higdon, J.J.L. The generation of feed-
The no-slip surface has a strong effect on the far-field behaviour of ing currents by flagellar motions. J.
Fluid Mech., 94:305–330, 1979
the flow. For example, for a vertical stokeslet, the flow falls off as
1/r3 , while u ∼ 1/r2 for a force parallel to the surface. The slower
fall-off for forces oriented parallel to the surface suggests, for exam-
ple, that in the action of cilia at the surface of organisms the move-
ments parallel to the surface are more important for propulsion and
fluid transport than those perpendicular to the surface.
The general result that gradients of elementary flow fields yield
higher-order singularities can be used to understand the flow field
generated by a bacterium swimming parallel to a no-slip surface.
Using the general formulation (6.146), we first obtain the flow field
in the x − y plane associated with a point force along x,
1 x2 x2 3x2 z
∥ F 1 z
ux = + 3 − + 3 − 2h + ,
8πµ r r R R R3 R5
(6.148a)
∥ F xy xy 6hxyz
uy = − 3+ . (6.148b)
8πµ r3 R R5
cellular motion (lectures 17-20) 137
The streamlines of this flow field are shown in Fig. 6.35(a), and we
observe from (6.148) that u ∼ 1/r2 in the far field. The gradient
of (6.148) is shown in Fig. 6.35(b), which is proportional to the flow
field of a stresslet near a no-slip surface. Again we see the generation
of strong recirculation regions.
of a stresslet in free space (6.99) to write the solution for the entire
flow field for a stresslet oriented along ex at position (0, 0, h) with
image at (0, 0, −h) near a free surface in the x − y plane as
2 2
P
∥ 3x 1 3x 1
u x ( x, z) = − x + − x , (6.149a)
8πµ r5 r3 R5 R3
2 2
P
∥ 3x 1 3x 1
uz ( x, z) = − ( z − h ) + − ( z + h ) ,
8πµ r5 r3 R5 R3
(6.149b)
p p
where r = x2 + (z − h)2 and R = x2 + (z + h)2 . It is straight-
∥ ∥
forward to verify that uz vanishes at z = 0 as does ∂u x /∂x there, as
required. If we now evaluate this flow at the stresslet location (0, h)
we find the wall-induced velocity is
∥ P
uz (0, h) = − . (6.150)
32πµh2
This is an attractive interaction for a pusher stresslet P > 0, and a
repulsion for a puller. A nearly identical result holds for a no-slip
surface, with only a change in the numerical factor in the denomi-
nator. Experiments31 indicate that this is likely the dominant effect 31
Berke, A.P. and Turner, L. and Berg,
leading to surface accumulation of bacteria and spermatozoa. H.C. and Lauga, E. Hydrodynamic At-
traction of Swimming Microorganisms
Next we consider a stresslet oriented at some angle ϕ with respect by Surfaces. Phys. Rev. Lett., 101:038102,
to the stress-free surface, for which the velocity components are 2008
dX 3 X 0
=− . (6.154)
dT π ( X 2 + 4)5/2
There are several interesting points about this dynamical law. First,
V (X)
it is a gradient flow in the sense that dX/dT = −dV ( X )/dx, with
1 1
V (X) = − , (6.155)
π ( x2 + 4)3/2 !1=8:
0 2 4 6
a potential illustrated in Fig. 6.38. We can therefore interpret the X
Figure 6.38: Effective potential for hy-
fluid-mediated interaction as a direct mechanical interaction between drodynamic interaction between two
the spheres, characterised by a potential energy function. stokeslets near a no-slip surface. Con-
tact at occurs at X = 2 (dashed line).
X
5
1
-4,000 -2,000 0
T
1 k 2 k
E+S −
↽−
−⇀− C −−→ E + P (7.1)
k−1
de dc
+ = 0, (7.4)
dt dt
and hence e + c = e0 at all times. The system of three remaining
reactions after eliminating that for p can thus be reduced to two by
setting e = e0 − c, yielding
ds
= − k 1 e0 s + ( k 1 s + k −1 ) c (7.5)
dt
dc
= k 1 e0 s − ( k 1 s + k −1 + k 2 ) c (7.6)
dt
As in all problems of this type, we now nondimensionalise the
equations. The enzyme is typically present at much lower concen-
trations than the substrate, so we identify the small parameter ϵ =
e0 /s0 . Let the rescaled time be τ = k1 e0 t, and set
s(t) k2
u(τ ) = , λ= , (7.7)
s0 k 1 s0
c(t) k + k2
v(τ ) = , k= 1 , (7.8)
e0 k 1 s0
Then we have the final equations for the substrate and complex,
du
= −u + (u + k − λ)v, u(0) = 1, (7.9a)
dτ
dv
ϵ = u − (u + k)v, v(0) = 0. (7.9b)
dτ
In this form, with ϵ multiplying the highest (time) derivative in
the problem, it is clear that a finite ϵ is a singular perturbation to the
kinetics & pattern formation (lectures 21-24) 143
0.5 Figure_701_and_704_MichaelisMenten.m.
v(= )(!)
v0 (= )(!!)
0
0 1 2 3
=
τ = (1 − u0 − k ln u0 ) /λ. (7.12) 1
R/Vmax
This satisfies the initial condition u0 (0) = 1, but the partner relation
(7.10) does not satisfy v0 (0) = 0. Nevertheless, the comparison in
Fig. 7.1 between u0 and v0 and the full numerical solutions shows
that they agree well beyond times of order ϵ. 1 s/Km
Returning to the fundamental result (7.10), we see in dimensional Figure 7.2: Michaelis-Menten reaction
rate versus substrate concentration.
units that the rate of reaction R = dp/dt would be
s
R = Vmax , (7.13)
s + Km
where Km = (k −1 + k2 )/k1 is known as the Michaelis constant and
Vmax = k2 e0 . As shown in Fig. 7.2, this rate law exhibits a linear
144 biological physics and fluid dynamics
u0 + k ln u0 = −λτ + C, (7.14)
where δ(ϵ) is the as yet unknown inner time that we will find on
the basis of balancing terms in the rescaled dynamics. By direct
substitution into (7.9a) and (7.9b) we obtain
1 dU
= −U + (U + k − λ)v, U (0) = 1, (7.16a)
δ dT
ϵ dV
= U − (U + k)V, V (0) = 0. (7.16b)
δ dT
where the inner solution is (U, V ). It is now clear that in order that
the time derivative balances the terms on the r.h.s. of (7.16b) we must
have δ(ϵ) ∝ ϵ and we are free to set the proportionality constant
to unity and thus find T = τ/ϵ. In other kinds of problems the
balancing of terms can be more subtle, but as ϵ only appears in the
time derivative, here it is simple. We therefore find that the equations
to be solved to find the inner solution are
dU
= ϵ [−U + (U + k − λ)v] , U (0) = 1, (7.17a)
dT
dV
= U − (U + k)V, V (0) = 0. (7.17b)
dT
kinetics & pattern formation (lectures 21-24) 145
dU0
= 0, U0 (0) = 1, (7.18a)
dT
dV0
= U0 − (U0 + k)V0 , V (0) = 0. (7.18b)
dT
This is easily solved to yield
1 h i
U0 ( T ) = 1, V0 ( T ) = 1 − e−(1+k)T . (7.19)
1+k
ψ(ϵ)
lim ψ(ϵ) = 0, lim = ∞. (7.20)
ϵ →0 ϵ →0 δ ( ϵ )
The matching condition is then simply that the limit of the outer
solution at some fixed Ti as ϵ → 0 equals that of the inner solution,
h i
lim u0 (τ )|τ =ψTi = lim U0 ( T )| T =Ti ψ/δ , (7.21)
ϵ →0 ϵ →0
0.5 Figure_701_and_704_MichaelisMenten.m.
vunif orm (= )
0
0 1 2 3
=
τ 1
ṽ0 (τ ) = v0 (τ ) + V0 − (7.24a)
ϵ 1+k
u0 ( τ ) e−(1+k)τ/ϵ
= − , (7.24b)
u0 ( τ ) + k 1+k
These results are shown in Fig. 7.4b, where we see excellent agree-
ment with the full numerics over the entire time range.
fractional binding
tional occupancy of binding sites on hemoglobin. As a function of
the partial pressure of oxygen in the blood, the occupancy exhibits a
characteristic "sigmoidal" shape shown in Fig. 7.5. This arises from
the feature of "cooperativity", in which the binding of one oxygen
pO2
molecule triggers a conformational change in the hemoglobin that
Figure 7.5: Schematic sigmoidal bind-
facilitates the binding of subsequent molecules. Such behaviour is ing curve of oxygen to hemoglobin, as
also known as allosterism and hemoglogin is said to be "allosteric". a function of the partial pressure of oxy-
This produces a very large change in absorbed oxygen over a rel- gen.
k1
−−
E+S ↽−⇀
− C1 (7.25a)
k−1
2k k3
C1 + S −−→ C2 −−→ C1 + P, (7.25b)
7.1.4 Slaving
The appearance of strong nonlinearities in the effective dynamics of
slow variables in systems with a strong separation of time scales can
be seen within a very simple model with only quadratic underly-
ing nonlinearities. Suppose that there are two degrees of freedom
p, q, where q is a slow variable (as in the Michaelis-Menten substrate
concentration) and p is fast (as in the MM complex concentration),
dq
= αq − βpq (7.29a)
dt
dp
ϵ = γp − δq2 . (7.29b)
dt
We interpret the terms αq and γp as representing autocatalysis. In
the steady state approximation, p = δq2 /γ and the slow dynamics is
dq βδ
≃ αq − q3 . (7.30)
dt γ
Importantly, the r.h.s. of (7.30) has the form −dV (q)/dq, where
1 βδ 4
V (q) = − αq2 + q . (7.31)
2 4γ
Thus, just as in our study of Euler buckling, we have a potential
capable of producing a pitchfork bifurcation as a function of α. This
corresponds to the physical scenario of overdamped dynamic, and in
fact produces a bifurcation as a function of α, from a function with
a single minimum for α < 0 to a double well. The key point to take
away from this simple example is that by slaving fast variables to
slow ones we can obtain very strong nonlinearities in the effective
equations of motion for the slow degrees of freedom, and these can
exhibit bifurcations and even bistability.
148 biological physics and fluid dynamics
ion channels
ion pumps
water + ions (low Na+, high K+, ...)
membrane
a few nanometers thick that separates the fluid inside from that sur-
rounding the axon. Within the membrane are specialized protein
structures that produce and regulate a large difference in concentra-
tions of certain ions in those fluid regions. These include "pumps"
that actively and selectively transfer ions across the membrane as
they consume chemical energy, and "channels" that can open and
close to regulate the passive flow of ions across the membrane. Na+
ions may have a concentration of 150 mM outside and 15 mM inside,
while K+ is 5 mM inside and 150 mM outside. As mentioned in the
introduction, the Nernst equation relates the electric potential differ-
ence across a membrane to the concentrations on the two sides
kB T Cout
V= ln (7.32)
e Cin
where the prefactor has the value 25 mV. The typical concentration
ratios yield voltages on the scale of 60 − 70 mV.
One of the key developments of HH, and Cole & Cole before
them,8 was to take advantage of the large size of the squid axon 8
K.S. Cole. Membranes, Ions and Im-
to insert an electrode inside it, oriented parallel to the long axis. This pulses. A Chapter in Classical Biophysics.
University of California Press, Berkeley,
allowed measurements of the electrical properties without any spa- 1968
tial effects (a "space clamp"), so that the unusual properties of the ion
channels could be studied in isolation. The basic phenomenology of
electrical impulse generation in an axon that was found is depicted
in Fig. 7.7. In the absence of any stimulus, there is a "resting poten-
tial" difference of ∼ 70 mV across the membrane. If a perturbation
is made, say by injecting current into the system, that voltage can
kinetics & pattern formation (lectures 21-24) 149
40
0 t (ms)
0 1 2 3 4 5
sation
repola
depolari
risatio
n
-55 threshold
-70 resting potential
refractory
period
hyperpolarisation
1
u̇ = cv, v̇ = − u, (7.38)
c
which is just the Hamiltonian dynamics u̇ = dH/dv, v̇ = −dH/du of
a harmonic oscillator with
1 2 1 2
H= u + cv . (7.39)
2c 2
Equivalently, ü = −u, so this structure naturally leads to oscilla-
tions at unit frequency. It is the competition between the dissipative
and oscillatory tendencies of the dynamics that leads to interesting
behaviour.
Finally, note that if we rescale time in (7.34) and (7.35) as t = cτ,
then they take the form
1
ϵu′ = v + u − u3 − V0 , (7.40a)
3
v′ = − (u − a + bv) , (7.40b)
v
0
These are plotted in Fig. 7.8. The slope of the cubic nullcline at
u = 0 is −1, whereas the second nullcline slope is −1/b, and thus -1
the restriction 0 < b < 1 guarantees the latter is larger and there is
-2 -1 0 1 2
only a single intersection of the two nullclines, and therefore only u
one fixed point (FP) (u∗ , v∗ ) and it can easily be shown that u∗ > 1. Figure 7.8: Nullclines of the FHN
If we write the generic two-component dynamics as model for V0 = 0. Obtained with
Figure_708_to_711_FitzHughNagumo.m.
where f u = ∂ f /∂u|u∗ ,v∗ that determines the stability of the FP. Per-
turbations of the form u = u∗ + δu, v = v∗ + δv will grow as eσt with
the growth rate σ found as the solution to the condition
fu − σ fv
= 0, (7.44)
gu gv − σ
1n p o
σ± = Tr ± Tr2 − 4Det . (7.45)
2
For stability, the real part of both roots must be negative, which re-
quires Tr < 0. If 0 < Det < Tr2 /4 the square root is real but less
then |Tr| and the roots are negative, while for larger Det the square
root yields an imaginary contribution, and still the real part of σ is
negative. Thus, stability requires Tr < 0 and Det > 0. In our case,
! !
c (1 − u ∗2 ) c − +
J = ∝ (7.46)
z =0 −1/c −b/c − −
2
(a) Figure 7.10: Dynamics of the FHN
2 (b) u
model II. As in Fig. 7.9 except
1
for V0 = 0.3, producing a sin-
1
gle firing event. Obtained with
Figure_708_to_711_FitzHughNagumo.m.
u; v
0
v
0
v
-1
-1
-2
-2 -1 0 1 2 0 10 20 30 40 50
u t
2
(a) Figure 7.11: Dynamics of the
2 (b) u
FHN model III. As in Fig. 7.9
1
except for V0 = 0.4, leading to
1
periodic firing. Obtained with
Figure_708_to_711_FitzHughNagumo.m.
u; v
0
v
0
v
-1
-1
-2
-2 -1 0 1 2 0 10 20 30 40 50
u t
We thus see that the FHN model in the slow inhibitor limit dis-
plays all the features of an "excitable medium" as exemplified by
neurons; small-amplitude response for stimuli below a threshold, a
single large-amplitude response beyond that threshold, and periodic
large-amplitude responses for even higher stimuli. An important les-
kinetics & pattern formation (lectures 21-24) 153
son from the study of this model is that the precise algebraic form
of the nullclines is less important than the overall geometry. The ex-
citable dynamics fundamentally requires nullcline intersection with
a single stable fixed point. In the next section we consider the situa-
tion with multiple fixed points and a fast inhibitor.
∂F
ut = mu xx + f (u) , where f (u) = − , (7.47)
∂u
F(u)
where we call the diffusion constant m for later convenience. By
construction, the dynamics of u is variational, with the form
Z ∞
δE
1 2
ut = − , E [u] = dx m (∇u) + F (u) . (7.48)
δu −∞ 2
In writing F this way, we have split out the first term, which repre-
sents the product of parabolas centred at u = 0 and u = 1, giving two
minimum with F (0; r ) = F (1; r ) = 0, and the second term that rep-
resents the contribution of the control parameter to changing those r 1 u
depths. The result is a function such that F (0; r ) = −(r − 1/2)/12
and F (1; r ) = (r − 1/2)/12 and the difference between the two is
1 1
∆F = F (1; r ) − F (0; r ) = r− . (7.51)
6 2
For r < 1/2 the state u = 1 is the more stable, and for r > 1/2 the Figure 7.14: The function f (u).
state u = 0 is more stable.
We first need to establish some terminology and methods that will
be useful for many subsequent discussions. The chief technique is a
linear stability analysis of a stationary state. Let us examine in turn the
three states ū that are zeros of f (u), and hence satisfy ut = 0. In each
case we write u = ū + û, linearize the dynamics in û, set û = eikx eσt
and find the growth rate σ (k) as a function of the wavevector k. If
σ < 0, ū is stable, while if σ > 0 ū is unstable. In adopting this
plane wave perturbation we are ignoring any lateral boundaries and
the boundary conditions (Dirichlet or Neumann) that might pertain
kinetics & pattern formation (lectures 21-24) 155
there. Such simplifications are not always justified (see Sec. 7.4.3),
but here they are helpfully simplifying. σ(k)
Near ū = 0, we find
ût = mû xx − (1 − r )û, (7.53a) Figure 7.15: Growth rate for the stable
state ū = 0.
2
σ = −(1 − r ) − mk , (7.53b)
u( x, t) = U ( x − vt), (7.59)
Figure 7.18: Generic form of F (u).
for some unknown speed v. Substituting, we find
where the energy difference ∆F is given in (7.51) and the first term
vanishes by the boundary conditions at infinity. We can then for-
mally solve for the front speed as
−∆F
v = R∞ 2
. (7.62)
−∞ dzUz
7.4.1 Phenomenology
Systems of reacting and diffusing chemicals are ubiquitous in bi-
ological systems, essentially responsible for life’s processes inside
cells. While many biological pathways involve numerous coupled
reaction-diffusion systems, much of our intuition is built on mod-
els with just a few coupled species. In this section we review on σ(k)
a phenomenological level some considerations that indicate the im-
portance of examining at least two coupled species in order to obtain
pattern-forming dynamics consistent with observations. This leads
k
naturally into Turing’s analysis, covered in the next section.
Let us start from a generic structure of a nonlinear PDE that em-
bodies diffusion and nonlinear chemical kinetics. From our previous
discussions we expect the form
A more interesting possibility is when both long and short wave- σ(k)
length are damped and σ(k ) is positive only in an interval of finite
k, as in Fig. 7.22. This band of unstable modes typically increases k
k*
in width as a control parameter is varied, starting from a critical k∗
when σ has a maximum just crossing the axis. The maximum of σ
defines the fastest growing mode, which often determines the wave-
length λ∗ = 2π/k∗ of the finite-amplitude patterns beyond onset.
In order to obtain growth curves like those in Fig. 7.22 with a sin- Figure 7.22: Growth curves for a
gle species we would require the linear operator L to have higher- pattern-forming system as a parameter
order derivatives than those from diffusion. Respecting parity sym- is increased (green arrow).
ut = αu − βu xx − γu4x + . . . . (7.66)
ut = D ∇2 u + f (u) − ρv (7.67a)
2
ϵvt = ∇ v + αu − βv, (7.67b)
in which we have scaled space and time in order that the diffusion
constant of the inhibitor is unity. The presence of the relative diffu-
sion constant D can lead to the key phenomenon of lateral inhibition
for D ≪ 1 as we discuss below. We may be interested in a whole
range of values for ϵ, not necessarily small. The various terms on
the right hand sides of these equations have the following interpre-
tation: The function f (u) embodies autocatalysis and bistability in
the manner we have discussed earlier in this chapter. The term −ρv
represents inhibition of the activator due to the presence of the in-
hibitor. In the inhibitor equation, αu represents stimulation of the
inhibitor due to the activator, while − βv describes self-limitation of
the inhibitor. Depending on the structure of f and the values of
the various coefficients, the FHN model can produce homogeneous
states, stripes and other periodic patterns, and even rotating spiral
waves. It is a "standard model" of pattern formation.
kinetics & pattern formation (lectures 21-24) 159
; = 0:2
-0.5 0.5 0.5 agram around u = 1. Obtained with
Turing Turing Figure_722_FHNgrowth.m.
stable stable
-1 1 1
0 5 10 0 0.5 1 0 0.5 1
k r r
Given the shape of σ (k) it is clear that the threshold for the onset of
a finite-wavelength instability occurs with σ (k) = 0 and dσ (k)/dk =
0 for some k c . The second relation may also be chosen as dσ/dk2 = 0
for convenience. One finds the threshold lies along a curve in the
r − ρ plane given by
√ √ 2
ρ c (r ) = r+ D , (7.77)
and
1 (black)
v0xx − v0 ∼ −u0 = − (7.82)
0 (white).
U (η ) ∼ U 0 (η ) + ρU 1 (η ) + · · · . (7.88)
j √
1
Qt ∼ ∓6 2D ∆F + ρ v0 ( x = Q j ) − . (7.93)
2
This general result can now be applied to several important cases.
For the case of a single front, we consider the case of a pattern that
is white (u = 0) to the right of an interface at x = Q(t) and black
(u = 1) to the left,
1, x ≤ Q(t)
u0 ( x ) ∼ (7.94)
0, x ≥ Q(t).
5
(b) Figure 7.24: Dynamics of two nearby
(a) fronts in the FHN model. Param-
4
eters are D = 0.04, r = 0.4 and
u(x; t); v(x; t)
Q(t)
hibitor (red) fields for two approach-
2 ing fronts, with time increasing up-
wards. Black circles indicate location
1 of u = 1/2, the midpoint of the fronts.
1 (b) Distance between fronts as a func-
0 0 tion of scaled time. Obtained with
-10 -5 0 5 10 0 5 10 Figure_724_FHNfronts.m.
p
x (1=2 ! r) 2Dt
Figure 7.24 shows the time evolution of a pair of fronts that start
at a large separation Q, with r < 1/2, so the u = 1 state is more
stable and the two fronts approach one another. As they do, the
leftward (negative) speed from the term r − 1/2 will be reduced as
Q decreases by the repulsive effect of the overlapping inhibitor fields,
eventually vanishing at a distance
1/2 − r
1
Q∗ = − ln . (7.101)
2 3ρ
pt = f u p + f v q + Du ∇2 p, (7.104a)
2
qt = gu p + gv q + Dv ∇ q, (7.104b)
While the trace of Jmod is just the sum of the traces of the two matri-
ces in (7.107),
Trmod = f u + gv −( Du + Dv )k2 , (7.108)
| {z }
Tr<0
and is < 0 since f u + gv < 0 by hypothesis and the diffusive contri-
butions are manifestly negative. In contrast, it is generally the case
for two matrices M and N that
Detmod = f u gv − f v gu − ( Du gv + Dv f u ) k2 + Du Dv k4 , (7.110)
| {z } | {z } | {z }
Det>0 ? >0
Du gv + Dv f u > 0, (7.111)
Det(k)
for then, as shown in Fig. 7.26, the function Detmod (k) can dip below
the axis. If we let x = k2 , then any crossings of the axis must be
roots of the quadratic equation Ax2 − Bx + C = 0, namely x± =
√
(1/2)( B ± B2 − 4AC ), with B2 − 4AC > 0. The smallest value of B
is that which yields the tangent condition shown as the middle curve
√
in Fig. 7.26, namely B = 4AC.
k
We thus conclude that the onset of the instability occurs at k*
p
Dv f u + Du gv = 4Du Dv Det. (7.112)
Figure 7.26: Behaviour of the determi-
nant of the stability matrix.
If we define the diffusivity ratio
Du
d= , (7.113)
Dv
then the condition (7.112) is
√
f u + dgv = 2 dDet. (7.114)
Since the r.h.s. of this equation is clearly > 0, we conclude that the
conditions that the homogeneous fixed point be stable and that the
diffusive system be unstable require the simultaneous equations
f u + gv < 0, (7.115a)
f u + dgv > 0, (7.115b)
Det 1/2
r
C
kc = = . (7.117)
A Du Dv
We may now summarize the fundamental features of the Turing
instability. A two-component reaction-diffusion system that has a
stable homogeneous FP can, given suitable conditions on the ratio of
diffusion constants of the two species, become unstable to a finite-
wavelength instability. As we have focused only on the results ob-
tainable from linear stability analysis we can say nothing definite
about the fate of a system beyond the instability, but generically one
expects the nonlinearities in the system to stabilize finite-amplitude
patterns with a characteristic wavelength close to that at onset. As
in many other pattern-forming systems,14 the precise form of those 14
Cross, M.C. and Hohenberg, P.C. Pat-
nonlinearities will determine the symmetries of the patterns, distin- tern formation outside of equilibrium.
Rev. Mod. Phys., 65:851–1112, 1993
guishing for example between stripes and spots.
biological physics). It was not until the early 1990s that experimen- 16
Castets, V. and Dulos, E. and Bois-
sonade, J. and De Kepper, P. Experi-
tal demonstrations of his predictions were reported. These were re- mental evidence of a sustained stand-
ported by the Bordeaux group of Castets et al.,16 , 17 followed by the ing Turing-type nonequilibrium chem-
Austin group of Ouyang and Swinney.18 Each of these groups uti- ical pattern. Phys. Rev. Lett., 64:2953–
2956, 1990
lized experimental setups with two reservoirs, each containing a non- 17
De Kepper, P. and Castets, V. and Du-
reacting subset of chemicals, linked by a porous medium in which los, E. and Boissonade, J. Turing-type
chemical patterns in the chlorite-iodide-
reactions occur without the disturbing effects of fluid motion. The
malonic acid reaction. Physica D, 49:
Austin experiments in particular allowed for spatially extended pat- 161–169, 1991
terns in two dimensions without any imposed chemical gradients. 18
Ouyang, Q. and Swinney, H.L. Tran-
sition from a uniform state to hexago-
Figure 7.27 shows that patterns in the form of stripes and hexagonal
nal and striped Turing patterns. Nature,
arrays of spots were observed. 352:610–612, 1991
We close this section by noting that reaction-diffusion patterns are 19
K.S. Kiprijanov. Chaos and beauty
not restricted to the stationary ones we have focused on. Turing in a beaker: The early historyof the
Belousov-Zhabotinsky reaction. Ann.
himself in his 1952 paper established that there could be oscillatory
Phys. (Berlin), 3-4:233–237, 2016
patterns as well, and we know this is the case from studies of the
Belousov-Zhabotinski reaction,19 which is famous for the rotating
spiral waves and target patterns that it produces (Fig. 7.28). Similar
patterns occur in a range of living systems, from slime molds to the
heart, as we mention below. The remarkable history of the discovery
of oscillatory dynamics in the BZ reaction also provides a case study
in the resistance of science to new ideas.20
(a) (b)
nt = −∇ · Jn , (7.118)
where the flux Jn of cells has both diffusive and chemotactic contri-
butions, written as
Jn = − Dn ∇n + nr ∇c. (7.119)
The first term is the standard Fickian form, where Dn is the cellu-
lar diffusion constant that arises from their random motions, and is
assumed constant. The chemotactic part embodies the fact that cells
climb up gradients in cAMP, and we can interpret the product r ∇c
as the typical chemotactic velocity uc , with the coefficient r (c) poten-
tially varying with the cAMP concentration. Keller and Segel suggest
a form r (c) ∼ 1/c that represents logarithmic sensing (uc ∼ ∇ ln c).
Thus, the variable n evolves according to the nonlinear equation
ct = Dc ∇2 c + n f (c) − k1 cp + k −1 i, (7.122a)
2
it = Di ∇ p + k1 cp − (k −1 + k2 ) i, (7.122b)
2
pt = D p ∇ p + ng(c, p) − k1 cp + (k −1 + k2 ) i. (7.122c)
k1 cp − (k −1 + k2 ) i = 0, (7.123)
removes two equations from the problem, leaving only the coupled
dynamics of the cell and cAMP concentrations, in the form
nt = Dn ∇2 n − ∇ · (rn∇c) , (7.124a)
2
ct = Dc ∇ c + n f (c) − k(c)c, (7.124b)
where
p0 k 2
k(c) = , (7.125)
Km + c
and Km = (k −1 + k2 )/k1 as in Michaelis-Menten theory.
There are obvious similarities between this two-variable system
and that of the Turing instability in the sense that the reaction terms
on the r.h.s. of (7.124b) embody production and degradation in fa-
miliar forms. But the key difference is the presence of an off-diagonal
spatial term from chemotaxis. To explore the consequences of this
term, we simplify the model further by taking the functions k (c) and
f (c) to be constant, and also assuming the chemotactic response co-
efficient r to be independent of c, yielding
n t = Dn ∇ 2 n − r ∇ · ( n ∇ c ) , (7.126a)
ct = Dc ∇2 c + αn − βc, (7.126b)
Noting that the n dynamics is of a flux form that conserved the total
number of cells, we choose the mean cell concentration n0 as a ref-
erence value for rescaling. The homogeneous system has an obvious
fixed point at
α
c0 = n0 . (7.127)
β
Setting r
n c α
ρ= , χ= , T = αt, X= x, (7.128)
n0 c0 Dc
we obtain
ρ T = d ∇2 ρ − γ ∇ · ( ρ ∇ χ ) , (7.129a)
2
χ T = ∇ χ + λ (ρ − χ) , (7.129b)
where
Dn rc0
d= , γ= . (7.130)
Dc Dc
With this simplified system we linearise around the fixed point
ρ = χ = 1 via
ρ = 1 + δρ, χ = 1 + δχ, (7.131)
to obtain the linearised system
−dk2 − σ γk2
2 = 0. (7.133)
λ −k − λ − σ
172 biological physics and fluid dynamics
1
σ+ (k) = (γ − d) k2 − 1 − d + 2γ2 k4 + · · · . (7.136)
λ
The same calculation shows that σ− is strictly negative as k → 0.
Thus, we see that as γ increases beyond d the first instability is at
k = 0, but there is a peak in the growth rate for γ > d, at k∗ ∼
p
γ − d. Figure 7.32 shows σ+ (k) as γ increases.
0
<(k)
-0.1
-0.2
0 1 2
k
This instability represents the transverse instability seen in the
Dictyostelium system (Fig. 7.31). The mechanism behind the instabil-
ity is the essential feedback between cAMP production and chemo-
taxis; a local fluctuation in cell density n leads to a local increase in
cAMP production, producing a gradient in cAMP the points toward
the peak in n, followed by chemotaxis toward the peak, further in-
creasing n. The potential runaway effect is tamed by both cell and
cAMP diffusion, along with cAMP degredation. Under certain cir- 23
Childress, S. and Percus, J.K. Non-
cumstances, in certain spatial dimensions, a true singularity does linear Aspects of Chemotaxis. Math.
Biosci., 56:217–237, 1981
arise from the KS equations.23
8
Bibliography
van der Waals, J.D. Over de Continuiteit van den Gas- en Vloeistoftoes-
tand (On the continuity of the gas and liquid state). PhD thesis, Leiden
University, 1873.
Holstein, B.R. The van der Waals interaction. Am. J. Phys., 69:441,
2001.
Bussi, Y., Shimoni, E., Weiner, A., Kapon, R., Charuvi, D., Nevo, R.,
Efrati, E., Reich, Z. Fundamental helical geometry consolidates
the plant photosynthetic membrane. Proc. Natl. Acad. Sci. USA,
116:22366–22375, 2019.
Sculley, M.J., Duniec, J.T., Thorne, S.W., Chow, W.S. and Boardman,
N.K. The stacking of chloroplast thylakoids. Quantitative analysis
of the balance of forces between thylakoid membranes of chloro-
plasts, and the role of divalent ions. Arch. Bioch. Biophys., 201:
339–346, 1980.
bibliography 175
Lis, L.J., McAlister, M., Fuller, N., Rand, R.P. and Parsegian, V.A. In-
teractions between neutral phospholipid bilayer membranes. Bio-
phys. J., 37:657, 1982.
Neuman, K.C. and Block, S.M. Optical trapping. Rev. Sci. Instrum.,
75:2787–2809, 2004.
Keller, J.B. and Rubinow, S.I. Slender-body theory for slow viscous
flow. J. Fluid Mech., 75:705–714, 1976.
Wiggins, C.H., Riveline, D., Ott, A. and Goldstein, R.E. Trapping and
wiggling: elastohydrodynamics of driven microfilaments. Biophys.
J., 74:1043–1060, 1998.
Guasto, J.S. and Johnson, K.A. and Gollub, J.P. Oscillatory Flows
Induced by Microorganisms Swimming in Two Dimensions. Phys.
Rev. Lett., 105:168102, 2010.
Leptos, K.C. and Chioccioli, M. and Furlan, S. and Pesci, A.I. and
Goldstein, R.E. Phototaxis of Chlamydomonas arises from a tuned
adaptive photoresponse shared with multicellular Volvocine green
algae. Phys. Rev. E, 107:014404, 2023.
Polin, M., Tuval, I., Drescher, K., Gollub, J.P. and Goldstein, R.E.
Chlamydomonas swims with two ‘gears’ in a eukaryotic version of
run-and-tumble locomotion. Science, 325:487–490, 2009.
Drescher, K., Leptos, K.C., Tuval, I., Ishikawa, T., Pedley, T.J. and
Goldstein, R.E. Dancing Volvox: Hydrodynamic Bound States of
Swimming Algae. Phys. Rev. Lett., 102:168101, 2009.
Berke, A.P. and Turner, L. and Berg, H.C. and Lauga, E. Hydrody-
namic Attraction of Swimming Microorganisms by Surfaces. Phys.
Rev. Lett., 101:038102, 2008.
Keller, E.F. and Segel, L.A. Model for chemotaxis. J. Theor. Biol., 30:
225–234, 1971.