Biofisica Teorica Testo Quantitativo Fondamentale PDF
Biofisica Teorica Testo Quantitativo Fondamentale PDF
Biofisica Teorica Testo Quantitativo Fondamentale PDF
ts by Ulrich
scrip Sc
re
tu
hw
Lec
arz
He it
y
id e rs
lb e rg U nive
2
Contents
Important numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Some history . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3
3.3 Charged wall in different limits . . . . . . . . . . . . . . . . . . . . 49
3.4 Poisson-Boltzmann theory . . . . . . . . . . . . . . . . . . . . . . . 52
3.5 Debye-Hückel theory . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.6 Strong coupling limit . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.7 Two charged walls . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.7.1 Poisson-Boltzmann solution . . . . . . . . . . . . . . . . . . 55
3.7.2 Debye-Hückel solution . . . . . . . . . . . . . . . . . . . . . 57
3.7.3 Strong coupling limit . . . . . . . . . . . . . . . . . . . . . . 58
3.8 Electrostatistics of viruses . . . . . . . . . . . . . . . . . . . . . . . 59
3.8.1 The line charge density of DNA . . . . . . . . . . . . . . . . 59
3.8.2 DNA packing in φ29 bacteriophage . . . . . . . . . . . . . . 60
3.8.3 Electrostatistics of viral capsid assembly . . . . . . . . . . . 63
4
5.4 Interacting polymers . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.4.1 Self-avoidance and Flory theory . . . . . . . . . . . . . . . . 127
5.4.2 Semiflexible polymer networks . . . . . . . . . . . . . . . . 129
7 Diffusion 151
7.1 Life at low Reynolds-number . . . . . . . . . . . . . . . . . . . . . 151
7.2 Measuring the diffusion constant . . . . . . . . . . . . . . . . . . . 156
7.2.1 Single particle tracking . . . . . . . . . . . . . . . . . . . . 156
7.2.2 Fluorescence Recovery After Photo-bleaching (FRAP) . . . 156
7.2.3 Fluorescence Correlation Spectroscopy (FCS) . . . . . . . . 160
7.3 Diffusion to capture . . . . . . . . . . . . . . . . . . . . . . . . . . 163
5
8.7 Basics of non-linear dynamics . . . . . . . . . . . . . . . . . . . . . 193
8.7.1 One-dimensional systems . . . . . . . . . . . . . . . . . . . 195
8.7.2 Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . 197
8.7.2.1 Saddle-node bifurcation . . . . . . . . . . . . . . . 198
8.7.2.2 Transcritical bifurcation . . . . . . . . . . . . . . . 198
8.7.2.3 Supercritical pitchfork bifurcation . . . . . . . . . 198
8.7.2.4 Subcritical pitchfork bifurcation . . . . . . . . . . 199
8.7.3 Two-dimensional systems . . . . . . . . . . . . . . . . . . . 200
8.7.4 Stable limit cycles . . . . . . . . . . . . . . . . . . . . . . . 202
8.8 Biological examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 204
8.8.1 Revisiting bistability . . . . . . . . . . . . . . . . . . . . . . 204
8.8.2 Stability of an adhesion cluster . . . . . . . . . . . . . . . . 205
8.8.3 Genetic switch . . . . . . . . . . . . . . . . . . . . . . . . . 207
8.8.4 Glycolysis oscillator . . . . . . . . . . . . . . . . . . . . . . 209
8.9 Excitable membranes and action potentials . . . . . . . . . . . . . 212
8.9.1 Channels and pumps . . . . . . . . . . . . . . . . . . . . . . 212
8.9.2 Hodgkin-Huxley model of excitation . . . . . . . . . . . . . 213
8.9.3 The FitzHugh-Nagumo model . . . . . . . . . . . . . . . . . 217
8.9.4 The cable equation . . . . . . . . . . . . . . . . . . . . . . . 220
8.9.5 Neuronal dynamics and neural networks . . . . . . . . . . . 221
8.10 Reaction-diffusion systems and the Turing instability . . . . . . . . 223
6
Important numbers
7
Some history (NP = Nobel Prize)
8
1990 Seifert and Lipowsky paper on vesicle adhesion
1991 spontaneous curvature phase diagram of vesicles (Seifert et al.)
1994 book Statistical Thermodynamics of Surfaces, Interfaces, and Membranes
by Safran
1994 area difference elasticity (ADE) model for vesicles (Miao et al.)
1995 Marko and Siggia model for stretching the WLC
1997 NP physics 1997 for laser cooling includes Steven Chu, who also works
on biomolecules
1997 RMP review by Jülicher, Armand and Prost on molecular motors
1998 MacKinnon Science paper on the structure of the K + channel (NP 2003)
2002 Lim et al. PNAS paper on shape of red blood cells
2014 NP chemistry for super-resolution microscopy to Eric Betzig, Stefan Hell
and Bill Moerner
2016 NP chemistry for the synthetic molecular motors (still missing is one on
biological molecular motors)
2018 NP physics for optical tweezers to Arthur Ashkin
9
10
Chapter 1
In this script on theoretical biophysics we will make use of concepts and methods
from many different fields of physics, which we will introduce when they are
needed. However, there are two parts of basic physics which we will need right
from the start, and therefore we briefly review them in this chapter. The first
one is statistical mechanics, and the second one is electrostatics.
In general, biological systems are not in equilibrium and driven by energy that
is supplied by the environment (food, light, etc). However, often state variables
change only slowly and therefore the system can be described by the laws of
thermodynamics and statistical physics, albeit often only on local and temporary
scales. Biological systems operate at relatively high and constant (body or room)
temperature and therefore the canonical ensemble is relevant.
We start by deriving the Boltzmann distribution for the canonical ensemble from
information theory. This approach to statistical mechanics has been pioneered
by Claude Shannon (founder of information theory) and Edwin Jaynes (inventor
of the maximum entropy principle). We start from the Shannon entropy
X
S=− pi ln pi (1.1)
i
where i numbers all states of the system and pi is the probability of a state with
P
i pi = 1. By multiplying with kB , we would get the physical (or Gibbs) entropy
S. For the microcanonical ensemble, we would have pi = 1/Ω being constant
(Ω is the number of states) and thus S = ln Ω, the famous Boltzmann formula.
One can show that entropy S is a unique measure for the disorder or information
content in the probability distribution {pi }.
We now want to maximize entropy under the constraint of constant average
P
energy, U = i Ei pi . We add normalization and average energy condition with
11
Lagrange multipliers:
X
δS = − (ln pi + 1 + α + βEi ) δpi = 0 (1.2)
i
leading to
pi = e−(1+α+βEi ) (1.3)
From the normalization we get
1
e−(1+α) = const = (1.4)
Z
with X
Z= e−βEi (1.5)
i
P
From the average condition U = (1/Z) i Ei e−βEi we get that β should be a
function of U . We can make the connection to temperature T and identify β =
kB T . Now we have the Boltzmann distribution:
1 −βEi
pi = e (1.6)
Z
where Z is the partition sum. For a continuous state space, we would replace
the sum over states by an integral over states. We conclude that the canonical
distribution is the one that maximizes disorder under the condition that the
average energy has a fixed (observed) value.
We now generalize to the case of particle exchange with a reservoir, for example
molecules in a bulk fluid that can adsorb or bind to a surface. Another example
might be the molecules in an open beer bottle lying on the floor of a lake. We
now have a second side constraint, namely for the average number of particles,
P
N = i Ni pi . Variation of the entropy gives
X
δS = − (ln pi + 1 + α + βEi + γNi ) δpi = 0 (1.7)
i
where we have introduced a third Lagrange parameter γ. With the same argu-
ments as above, we can identify γ = −βµ with the chemical potential µ. We then
get X
ZG = e−β(Ei −µNi ) (1.8)
i
12
1.1.3 The harmonic system
We now consider a system with one harmonic degree of freedom at constant
temperature. This could be for example a particle in a one-dimensional laser
trap with a harmonic potential E = 21 kx2 , where k is the spring constant (trap
stiffness) and x is the one-dimensional state space coordinate (position). The
corresponding partition sum is
∞ ∞ 1
1 2πkB T
ˆ ˆ
2
Z= dx exp(−βE) = dx exp(−β kx2 ) = (1.10)
−∞ −∞ 2 k
2
where β = 1/(kB T ) and we have evaluated the Gaussian integral dxe−ax =
´
σx2 =< (x− < x >)2 >=< (x2 − 2x < x > + < x >2 ) > (1.14)
kB T
=< x2 > − < x >2 =< x2 >= (1.15)
k
The average energy is
1 −1 kB T
ˆ
< E >= dx E exp(−βE) = ∂β Z = −∂β ln Z = (1.16)
Z Z 2
This is the famous equipartition theorem: every harmonic degree of freedom
carries an energy of kB T /2. Here we have one degree of freedom, for a harmonic
oscillator it would be two (potential and kinetic energy) and for an ideal gas
with N particles it would be 3N (only kinetic energy, but N particles in three
dimensions). The specific heat is constant:
kB
cV = ∂T < E >= (1.17)
2
For the variance of the energy we find
2 1 2 1
σE =< E 2 > − < E >2 = ∂β Z − ( ∂β Z)2 (1.18)
Z Z
(k T )2
B
= ∂β2 ln Z = −∂β < E >= (1.19)
2
13
A B C
F F
Figure 1.1: (A) A dielectric bead is attracted to the center of the laser beam. The force
F is proportional to the distance from the center. This is the principle of the optical
tweezer as developed in the 1970s by Arthur Ashkin at Bell Labs (Nobel prize physics
2018). (B) The optical tweezer can be used e.g. to measure the force-velocity relation
of a molecular motor. Using a feedback system that keeps force F constant, one can
measure the corresponding velocity v of the motor. For calibration of trap stiffness k,
one uses the relation < x2 >= kB T /k for a harmonic system. (C) Force-velocity relation
for the molecular motor kinesin (from Mark J. Schnitzer, Koen Visscher and Steven M.
Block, Force production by single kinesin motors, Nature Cell Biology 2, 718 - 723, 2000).
The larger F , the small v: eventually the motor gets stalled (v = 0) at high forces. The
higher ATP-concentration, the faster the motor.
where H = i p~2i /2m is the ideal gas Hamiltonian (only kinetic energy), p~i and
P
14
where s
h2
λ= (1.23)
2πmkB T
is the thermal (de Broglie) wavelength. A typical value for an atom is 0.1
Angstrom and below this scale, quantum mechanics become relevant. The free
energy follows with the help of Stirling’s formula ln N ! ≈ N ln N − N for large N
as
!
zN V
F = −kB T ln Z = −kB T ln = −kB T N ln +1 (1.24)
N! λ3 N
The Euler fundamental form for the Helmholtz free energy F = F (N, V, T ) is
From the statistical mechanics result for the free energy F = F (N, V, T ), we can
thus now calculate the pressure p as
N
p = −∂V F = kB T ⇒ pV = N kB T (1.26)
V
The result is known as the thermal equation of state or simply as the ideal gas
law.
The average energy is the caloric equation of state:
3N
< E >= −∂β ln Z = −N ∂β ln β −3/2 = kB T (1.27)
2
with p0 = kB T /λ3 (note that from the three terms, two have canceled each other).
Thus chemical potential and pressure are related logarithmically.
We can write our fundamental equation F (T, V, N ) and the three equations of
state in a very compact way using density ρ = N/V :
F
f= = kB T ρ ln(ρλ3 ) − 1 (1.29)
V
p = ρkB T (1.30)
<E> 3
e= = ρkB T (1.31)
V 2
µ = kB T ln ρλ3 (1.32)
15
1.1.5 The law of mass action and ATP-hydrolysis
From the ideal gas, we get for the chemical potential of species i in dilute solution:
ci
µi = µi0 + kB T ln (1.33)
ci0
Thus the change in Gibbs free energy at constant T and constant p is
X ∂G X X
∆G = ∆Ni = µi ∆Ni = µi νi ∆N (1.34)
i
∂Ni i i
where νi are the stoichiometric coefficients of the reaction and ∆N is the reaction
coordinate. At equilibrium, ∆G = 0 and ∆N drops out:
X ci,eq
0= νi µi0 + kB T ln (1.35)
i
ci0
This leads to
∆G = ∆G0 + kB T ln (Πcνi i ) , ∆G0 = −kB T ln Keq (1.38)
with the understanding that to get a dimensionless argument of the logarithm,
we might have to insert some reference concentration (typically 1 M). A very
important example is ATP-hydrolysis, for which we have νAT P = −1, νADP = +1
and νPi = +1. Thus we get
[ADP ][Pi ]
∆G = ∆G0 + kB T ln (1.39)
[AT P ]
With a reference concentration of 1M , the first term is −12.5kB T . For cellular
concentrations ([ADP ] = 10µM, [Pi ] = mM, [AT P ] = mM ), the second term is
−11.5kB T , so together we have ∆G = −24kB T .
16
where
VN
Zideal = (1.41)
N !λ3N
as above and
N ˆ
1 Y
Zinter = d~qi e−βU ({~qi }) (1.42)
V N i=1
This term does not factor into single particle functions because the potential U
couples all coordinates. Yet all thermodynamic quantities separate into an ideal
gas part and a correction due to the interactions. In particular, we have
F = −kB T ln Z = Fideal + Finter (1.43)
p = −∂V F = pideal + pinter (1.44)
The formulae for the ideal expressions have been given above. For the pressure,
one expects that the correction terms should scale at least in second order in
ρ, because two particles have to meet in order to give a contribution to this
term. This suggests to make the following ansatz of a Taylor expansion in ρ, the
so-called virial expansion:
∞
X
pinter = kB T Bi (T )ρi (1.45)
i=2
where the Bi (T ) are called virial coefficients. For many purposes, it is sufficient to
consider only the first term in this expansion, that is the second virial coefficient
B2 (T ). We then have
h i
F = N kB T ln(ρλ3 ) − 1 + B2 ρ (1.46)
p = ρkB T [1 + B2 ρ] (1.47)
For pairwise additive potentials, one can show
1
ˆ
B2 (T ) = − d~r e−βu(~r) − 1 (1.48)
2
For the van der Waals model, one considers two effect: a hard core repulsion with
particle diameter d and a square well attractive potential with an interaction
range δ and a depth ǫ. Then one gets
2π 3 ǫ a
B2 (T ) = d − 2π(d2 δ) =b− (1.49)
3 kB T kB T
where we have introduced two positive constants b (four times the repulsive eigen-
volume) and a (representing the attractive part). This general form of B2 (T ) has
been confirmed experimentally for many real gases. It now allows to rewrite the
gas law in the following way:
N
pV = N kB T (1 + B2 ) (1.50)
V
N N 2a
= N kB T (1 + b ) − (1.51)
V V
N kB T N a2
≈ N
− (1.52)
1 − bV V
17
a) b)
T T
F S F S
Tc
G L
Tt
ρ ρ
c) d)
ρ p
S
L F
Tt Tc T Tt Tc T
Figure 1.2: Combining the fluid-fluid and the fluid-solid phase transitions in (a), we get
the complete phase diagram of a simple one-component system in (b). In (c) we swap
T and ρ axes. By replacing ρ by p, we get the final phase diagram in (d). Two-phase
coexistence regions become lines in this representation. Carbon dioxide has a phase
diagram like this, but water does not.
thus
kB T a
p= − 2 (1.53)
(v − b) v
where v = V /N = 1/ρ is the volume per particle. This is the van der Waals
equation of state: the volume per particle is reduced from v to v − b due to
excluded volume, and pressure is reduced by the attractive interaction, that is
less momentum is transfered onto the walls due to the cohesive energy.
The van der Waals equation of state (1.53) is characterized by an instability. For
a stable system, if a fluctuation occurs to higher density (smaller volume), then
a larger pressure should result, which can counteract the fluctuation. Therefore
thermodynamic stability requires
∂p
<0 (1.54)
∂V
However, below the critical temperature Tc = (8a)/(27bk) the van der Waals
isotherms from (1.53) have sections in which this stability criterion is violated.
This indicates a fluid-fluid phase transition. The transition region can be cal-
culated by the Maxwell construction from thermodynamics. The van der Waals
gas thus predicts the fluid-fluid (gas-liquid) phase coexistence observed at low
temperatures.
18
Interacting systems also show a phase transition to a solid at high densities.
Together, one gets the generic phase diagram for a one-component fluid. It fits
nicely to the experimental results for simple fluids such as carbon dioxide (CO2 ).
However, the phase diagram for water is different, as we will see later.
The physics of phase transitions is important to understand how water molecules,
lipids and proteins behave as collectives. Traditionally, it has been thought that
biological systems tend to avoid phase transitions and that their appearance is
a sign of disease (like gall and bile stones). In recent years, it has become clear
however that biology uses phase transitions to create membrane-free compart-
ments, e.g. P-granules1 . A related and very important subject are rafts, which
are phase-separated but kinetically trapped domains in biological membranes2 .
1.2 Electrostatics
1.2.1 Electrostatic potential
1 ~r − ~r′
ˆ
~ =
E d~r′ ρ(~r′) ~
= −∇Φ (1.55)
4πǫ0 ǫ |~r − ~r′|3
1 ´ ρ(~r′)
Φ(~r) = d~r′ (1.56)
4πǫ0 ǫ |~r − ~r′|
~ ×E
∇ ~ =0 (1.57)
The Poisson equation implies that charges are the sources for the electrostatic
potential. For instance, the potential of a point charge with volume charge dis-
1
For a recent review, compare Hyman, Anthony A., Christoph A. Weber, and Frank Juelicher,
"Liquid-liquid phase separation in biology", Annual review of cell and developmental biology 30
(2014): 39-58.
2
For a review, compare Lingwood, Daniel, and Kai Simons, "Lipid rafts as a membrane-
organizing principle", Science 327.5961 (2010): 46-50.
19
tribution ρ(~r) = Q · δ(~r) can directly be calculated from equation 1.58:
spherical
symmetry 1 d2 (rΦ)
∇2 Φ = =0
r dr2
d(rΦ) A2
⇒ = A1 ⇒ Φ = A1 +
dr r
As an appropriate boundary condition we choose Φ(∞) = 0, hence A1 = 0. By
comparing our result with the Poisson equation, we finally get
Q 1
⇒ Φ(r) = ·
4πǫ0 ǫ r
which is the well-known Coulomb law. From a mathematical point of view, the
Coulomb law is the Green’s function (or propagator) for the Laplace
operator in
1
3D. The given solution can be checked to be true because ∆ r = −4πδ(r).
Sometimes it is useful to rewrite equation 1.58 in an integral form, using the
divergence theorem known from vector calculus:
divergence Poisson
theorem equation ρ(~r)
ˆ ˆ ˆ
~ ~
E dA = ~ ~
d~r ∇ · E = d~r
∂V V ǫ0 ǫ
~ = Qv
~ dA (1.59)
´
⇒ ∂V E Gauss law
ǫ0 ǫ
where ∂V is a closed surface, V its enclosed volume and QV the enclosed charge.
As an example equation 1.59 was used to compute the radial component Er of the
electric field of rotationally symmetric charge distributions (note that the angular
components vanish due to spatial symmetry). For a large sphere the Gauss law
reads:
~ = 4πr2 Er = QV ⇒ Er = QV
ˆ
~ dA
E
∂V ǫ0 ǫ 4πǫ0 ǫr2
thus we again recover Coulomb’s law, as we should.
The reference position ~r1 can be taken to be at infinity, where the potential
vanishes. For a continuous charge distribution, we therefore have
ˆ
Epot = d~r′ ρ(~r′ )Φ(~r′ )
20
We now consider a charge distribution localized around the position ~r and perform
a Taylor expansion around this point:
ˆ h i
~ r′ ) + . . . = QΦ(~r) − p~ · E
Epot = d~r′ ρ(~r′ ) Φ(~r) + (~r′ − ~r)∇Φ(~ ~ + ...
where the monopole Q is the overall charge and the dipole is defined as
ˆ
p~ = d~r′ ρ(~r′ )(~r′ − ~r)
We now write the interaction potential between two charge distributions. For a
monopole Q1 at the origin interacting with a monopole Q2 at ~r, we simply get
back Coulomb’s law:
Q1 Q2 1
Epot =
4πǫ0 ǫ r
by using the first term and the potential from a monopole. For a dipole p~ at ~r
interating with a monopole Q at the origin, we use the second term:
~ = Q p~ · ~r
Epot = −~
p·E
4πǫ0 ǫ r3
For two dipoles interacting with each other, we first take the potential resulting
from a dipole at the origin, which can be read off from the preceding equation:
1 p~1 · ~r
Φ=
4πǫ0 ǫ r3
We then get for the interaction
1 p~1 · p~2 3(p~1 · ~r)(p~2 · ~r)
Epot p2 · E~1 =
= −~ −
4πǫ0 ǫ r3 r5
The dipolar interaction is very prominent in biological systems. In particular,
water carries a permanent dipole and thus water molecules interact with this
potential function.
21
22
Chapter 2
From general physics, we know that there are four fundamental forces, but of
these only the electrostatic force is relevant for cohesion in biological material.
However, it comes in many different forms, as we shall see in this section. We start
with a discussion of the mechanical properties of biomaterial and immediately see
that we are dealing with very weak interactions on the order of thermal energy
kB T . We then review the details of these interactions and how they can be used
in molecular and Brownian dynamics simulations to predict the behaviour of
biomolecules, most prominently of proteins.
∆L
effect: strain ǫ= L [ǫ] = 1
23
Figure 2.1: Different ways to measure the mechanical rigidity of single cells. (a) Cell
stretching between two microplates. A related setup is pulling with an atomic force mi-
croscope (AFM). (b) Cell stretching with the optical stretcher. The functional principle
is similar to that of optical tweezers.
The simplest possible relation between the two quantities is a linear one:
σ =E·ǫ (2.1)
where E is the Young’s modulus or rigidity of the material with [E] = P a.
For cells, this elastic constant is typically in the order of 10 kP a. This is also
the typical stiffness of connective tissue, including our skin. In general, tissue
stiffness is in this range (on the cellular scale, the softest tissue is brain with 100
Pa, and the stiffest tissue is bone with 50 kPa).
Equation 2.1 might be recognized as Hooke’s law, and in fact we can think of the
macroscopic deformation as the effect of the stretching of a huge set of microscopic
springs which correspond to the elastic elements within the material. Equation
2.1 can be rewritten as
E·A
F = · ∆L (2.2)
L
thus k = E · A/L is the "spring constant" of the material. EA is often called the
1D modulus of the material.
Let us now assume that the system is characterized by one typical energy U and
one typical length a. A dimensional analysis of E yields E ≈ aU3 . As an example
a crosslinked polymer gel as illustrated in figure 2.2b can be considered.
The elasticity of cellular material is determined by supramolecular complexes
forming the structural elements of the cell with a typical scale a = 10 nm. There-
fore we get for the typical energy
U = E · a3 = 10kP a · (10 nm)3 = 10−20 J (2.3)
This is in the order of the thermal energy at ambient or body temperature (300 K)
known from statistical mechanics:
J
kB T = 1.38 · 10−23 · 300 K = 4.1 · 10−21 J = 4.1 pN nm (2.4)
K
J
where kB = 1.38 · 10−23 K is the Boltzmann constant.
In physical chemistry, one usually refers to moles rather than to single molecules:
kJ kcal
kB T · NA = R · T = 2.5 = 0.6 (2.5)
mol mol
24
Figure 2.2: (a) A slab of elastic material of length L and cross sectional area A is
stretched by a force F . The force acting on the material will result in a deformation. In
the case shown here the box will be stretched by the length ∆L. (b) Illustration of a
polymer gel with a meshsize a corresponding to a typical length scale. In this example,
the typical energy U is the elastic energy stored in a unit cell.
σ = η · ǫ̇ (2.6)
and a typical value for the viscosity of cells is η is 105 P a s, which is 8 orders
of magnitude larger than for water. This high viscosity comes from the polymer
networks inside the cell. The corresponding time scale is
and corresponds to the time the system needs to relax from the external per-
turbations by internal rearrangements. However, these consideration are only
relevant on cellular scales. If we make rheological experiments on the scale of
molecules, then we are back to the viscosity and relaxation times of water.
25
Chemical Bond Bond Energy
C −C 140 kB T
C=C 240 kB T
C≡C 330 kB T
H − CHO 144 kB T
H − CN 200 kB T
Table 2.1: Some chemical bonds and their corresponding bond energies (at T ≈ 300 K)
dU q1 q2
F =− ∼+ 2 (2.9)
dr r
which is repulsive if the two electric charges q1 and q2 have the same sign and
attractive otherwise.
The Coulomb interaction is a "long-ranged" interaction in 3D. To illustrate this,
consider the cohesive energy density of a bulk material of diameter L:
" #
L 3−n
ˆ L
2 1 3−n L 3−n
Utot = dr r n ∼ r |a = a −1 (2.10)
a r a
where a is a microscopic cutoff due to the Born repulsion. Taking the limit
L → ∞ in equation 2.10 shows that Utot does not diverge for n > 3, correspond-
ing to a short-ranged interaction where only the local environment significantly
contributes to the force on a point-like object. On the other hand, for n < 3 the
26
Figure 2.3: (a) Diameter of a typical cell in comparison to the thickness of the biological
lipid bilayer membrane. Note the very strong separation of length scales: a very thin
oily layer holds together the very large cell. (b) Drop of ǫ across the membrane. The
situation is similar to two metal sheets separated by plastic. Thus the membrane forms
a capacitor.
27
Figure 2.4: (a) Schematic drawing of a simple ionic crystal (such as NaCl). (b) Two
dipoles with dipole moments p~1 = e · a · n~1 and p~2 = e · a · n~2 , respectively.
For the total energy density of a crystal, one has to sum over all interactions
between nearest, next-nearest,... neighbours within the crystal. Let us first
consider only one row (compare figure 2.4a).
e2 2e2
1 1
Urow = · 2 · −1 + − + . . . =− ln 2 (2.11)
4πǫ0 a 2 3 4πǫ0 a
e2 N kcal
Utot = − 1.747 = −206 (2.12)
| {z } 4πǫ0 a mol
Madelung
constant
From the negative sign of the total energy in equation 2.12 it can be concluded
that the crystal is stable. The vaporization energy of a NaCl crystal was experi-
mentally determined to be 183 kcal
mol . Hence, although equation 2.12 is the result
of strong assumptions, it nevertheless agrees relatively well with the experimental
value.
(ea)2 h
ˆ n~2 · ~rˆ
i
U= n
~ 1 · n
~ 2 − 3 n
~ 1 · ~
r (2.13)
ǫ0 ǫr r3 | {z }
f (Θ,φ,... )
The factor f (Θ, φ, . . . ) does not depend on distance, but on all angles involved.
It is thus determined by the geometrical arrangement of the two dipoles and
its sign determines whether a certain orientation is favourable or not. Figure
28
Figure 2.5: f-values of different geometrical arrangements of two dipoles. The more
negative the f-value becomes, the more favourable is the arrangement.
2.5 shows some dipole arrangements and their corresponding f-values. The most
favorable orientation is a head-tail-alignment. In water, but also in dipolar fluids
and ferrofluids, this leads to dipole chains, network formation and spontaneous
polarization.
The interaction between charge distributions is further weakened by thermal mo-
tion. If the dipoles are free to rotate, the interaction becomes weaker. For
example, if a charge Q is separated by a distance r from a dipole with dipole
moment µ ~ = q · a, as depicted in figure 2.6, the electrostatic energy of the system
is given by
Qµ
U (~r, Θ) = · cos(Θ) (2.14)
4πǫ0 ǫr2 | {z }
| {z } orientation factor
U0
The dipole is rotating due to thermal forces, that is why we calculate an effective
interaction law by a thermal average weighted with the Boltzmann factor:
sin(Θ)dΘU (~r, Θ) exp −UkB(~rT,Θ)
´π
0
U (~r) = ´π
−U (~
r,Θ)
(2.15)
0 sin(Θ)dΘ exp kB T
−U (~
r,Θ)
If we assume that the interaction is weak compared to thermal energy, kB T ≪
1, then we can simplify the above expression:
−d(cos(Θ))U0 cos(Θ) 1 − U0kcos(Θ)
´π
0 BT
U (~r) = ´π
U0 cos(Θ)
0 −d(cos(Θ)) 1 − kB T
2
U02
1 Qµ 1
= − =− · (2.16)
3kB T 3kB T 4πǫ0 ǫ r4
So we see the change in the interaction potential from r12 for a static dipole to
1
r4
for a rotating one. The thermal motion weakens the Coulomb interaction also
for dipole-dipole interaction. A similar calculation can be made for dipoles that
are free to rotate with a centre-to-centre separation of r. We then obtain
2
2 µ1 µ2 1
U (~r) = − · (2.17)
3kB T 4πǫ0 ǫ r6
29
Figure 2.6: Interaction between a single charge and a rotating dipole.
30
O
H H
O
H H
Thus two permanent dipoles interact with an attractive and short-ranged 1/r6 -
potential.
A universal and short-ranged 1/r6 -attraction also arises for completely neutral
atoms due to quantum fluctuations. A neutral atom can always form a dipole by
quantum fluctuations, and this induces another dipole in a near-by atom, with
an interaction potential
31
Figure 2.9: Tetrahedral structure of ice and water, due to the hydrogen bonds between
the water molecules. Source: www.lsbu.ac.uk.
While the van der Waals interaction tends to condense water molecules, the
network of hydrogen bonds creates a more open structure. Because the second
effect dominates in ice, it swims on water. This also leads to the maximal density
of water at 4 C◦ . Pressure squeezes the molecules together and usually leads to
freezing; in water, it leads to melting. This is part of the explanation why you
can skate on ice, but not on glass. The feature of water is demonstrated in figure
2.10, where the phase diagrams of water and carbon dioxide are compared.
Water is also a very special solvent. It is ordered by the presence of the solutes.
For a hydrophobic solute, water molecules point their H-bonds away from the
solute. This decreases the entropy and therefore makes solution unfavorable
(measured by calorimetry, the effect is the strongest at 25 C◦ ). Because of the
“hydrophobic effect” water and oil do not mix. Non-polar solutes attract each
other in water and this phenomenon is called the “hydrophobic interaction”.
The large energy stored in the network of hydrogen bonds results in large values
for the surface tension, melting and boiling temperatures, heat capacity, etc.
Because the network of hydrogen bonds is easily polarized, water has a very
high dielectric constant (ǫ = 80). It is also important to remember that polar
solutes prefer polar solvents due to the low self-energy. In analogy to the previous
paragraph this effect is called “hydrophilic interaction”.
32
Figure 2.10: Phase diagrams of water on the left and of carbon dioxide on the right.
Source: www.astrosociety.org.
33
Figure 2.12: HP-model on an 3 × 2 - lattice. The upper panel shows the possible config-
urations on this lattice. Center panel: Possibilities to arrange the sequence HPHPHP on
the lattice. Note, that for the third configuration there exist two possible arrangements.
The energy penalty per H-P contact is ǫ (denoted as green lines). Recall that the envi-
ronment of the polymer is polar. Lower panel: PHPPHP sequence on the lattice. The
first configuration has a unique lowest energy and therefore forms the ground state.
The special properties of water are not only the basis of membrane assembly, but
also of protein folding. A standard method for studying and analyzing protein
folding is the HP-model by Ken Dill. It has been extensively studied on the
lattice by exact enumeration. The standard case is a 3 × 3 × 3 lattice, which can
contain 227 = 134217721 sequences and has 103346 possible configurations (this
number is non-trivial because one has to figure out all symmetry operations that
make two configurations identical in order to avoid overcounting). We pick one
configuration and fill it with a given sequence. After finishing the construct on
the lattice, for every amino acid positioned on the outside of the lattice an extra
P is added. After that every unfavorable H-P contact is assigned a free energy
penalty ǫ. This is repeated for all configurations and at the and we look for the
one with the lowest energy for a given sequence. If this ground state is unique,
we call it “native structure” and the sequence is “protein-like”.
The HP-model is a very useful toy model for protein folding. We now consider
a simplier variant. This time we have a 2 × 3 lattice with 26 sequences and 3
different configurations. The solvent molecules surrounding the lattice pattern
are assumed to be P-monomers. We now try to fit two different sequences on this
lattice — HPHPHP and PHPPHP. In figure 2.12 all possible configurations for
both sequences are shown.
While the first sequence (HPHPHP) is degenerated, the second (PHPPHP) has
a unique ground state. The sequence is therefore protein-like. The probability
to find the chain in the native structure as function of temperature is given by a
sigmoidal function, see figure 2.12:
exp(−2βǫ)
Pfold = (2.20)
exp(−2βǫ) + 2 exp(−4βǫ)
34
1.1
1
0.9
0.8
Pfold
0.7
0.6
0.5
0.4
0 1 2 3 4 5 6
TkB/ε
Figure 2.13: The probability to find the native structure as a function of temperature.
Figure 2.14: Polymer brushes as an example for steric interaction. As the brushes ap-
proach each other, the volume available for their motion and hence the entropy is reduced,
leading to an effective repulsion.
35
Figure 2.15: The planes on the left side represent two membranes that fluctuate to and
away from each other. On the right hand side one can see two whole cells with fluctuating
membranes .
The last example given here is the depletion interaction. Imagine two large
particles (depicted as large spheres in figure 2.16) that are surrounded by many
small particles. The volume available to the small molecules is marked blue, the
excluded volume is marked red. When the two large particles come close together,
so that the restricted volumes on their surfaces start to overlap, the entropy of the
system increases, because the volume available to the small molecules increases.
The system tries to reach a state with higher entropy, that is why the interaction
is called entropic attraction.
For a long time, it was thought that biological systems tend to avoid phase
separation, in the sense of the one-component phase diagram shown in Fig. 1.2,
which has two types of phase separations, fluid-crystal and fluid-fluid.
Crystal formation occurs in the human body, either as calcium phosphate crys-
tals in the bone, which we need, or as gall or bladder stones, which we do not
need. Bacteria or viruses sometimes use gene crystals to go into a dormant stage,
to survive for a long time under unfavorable conditions. However, all of these
examples are rather special and it is believed that usually concentration are not
sufficiently high to effect crystal formation. Obviously one artificially cranks up
36
concentration to get the protein cyrstals needed for structure determination by
X-ray diffraction.
Moreover, there was not much evidence for liquid-liquid phase separation (be-
cause the proteins are immersed in the cytoplasm, the fluid-fluid phase separa-
tion corresponds to low and high density solutions of proteins (and not to gas
and liquid), so we call them liquid-liquid in this context.) Recently, this notion
has changed completely, because it was realized that many proteins are disor-
dered rather than structured, and that intrinsically disordered proteins (IDPs)
tend to undergo liquid-liquid phase separations. Examples include the nucleolus,
P-granules and stress granules. It is now realized that liquid-liquid phase separa-
tion are a convenient way to establish compartments in the cell, as an alternative
to using membranes or protein shells to create closed compartments, like in vesi-
cles or viruses, respectively. With bacterial microcompartments (BMCs), there is
even an example which combines both, liquid-liquid phase separation and protein
capsid formation.
From the theoretical point of view, one can say that structured and disordered
proteins behave as colloids and polymers, respectively. Both phase separate, but
large colloids do only crystallize, while polymers only have a liquid-liquid loop.
Thus in each case, one of the two transitions from Fig. 1.2 is missing.
37
Figure 2.17: Energy distribution over time.
1 d2
T aylor
d
r~i (t)∆t2 ± ...
Expansion
r~i (t ± ∆t) = r~i (t) ± r~i (t) · ∆t +
dt 2 dt2
now we add both equations and get
F~i (t)
r~i (t + ∆t) = 2~
ri (t) − r~i (t − ∆t) + · ∆t2
mi
+O((∆t)4 ) (2.25)
One advantage is that the odd terms drop out, but more importantly the velocities
are not needed and can be calculated independently by
38
effects, one can work with periodic boundary conditions, truncated Lenard-Jones
potentials and the appropriate Ewald sum for the Coulomb interaction. The
ensemble described here is a NVE ensemble. If the temperature is fixed and the
energy fluctuates we have a canonical or NVT ensemble. The standard solution
in this case is the Nose-Hoover-thermostat. There are several classical papers
and book that describe the details of this important method.
Here is a list of the classical papers on this subject:
• Warshel, A., and M. Karplus. "Calculation of ground and excited state po-
tential surfaces of conjugated molecules. I. Formulation and parametriza-
tion." Journal of the American Chemical Society 94.16 (1972): 5612-5625.
• Car, Richard, and Mark Parrinello. "Unified approach for molecular dy-
namics and density-functional theory." Physical review letters 55.22 (1985):
2471.
In 2013, the Nobel prize for chemistry was awarded to Karplus, Levitt and
Warshel for the development of MD.
Books on MD-simulations:
39
• GROMOS: GROningen MOlecular Simulation computer program package
(Wilfred van Gunsteren, Switzerland)
mv̇ = F − ξv + ση(t)
Note that it is mandatory to add both terms together, because the noise term
alone would input too much energy into the system, so the damping is required
to balance this effect. This equation is the famous Langevin equation. It is a
stochastic differential equation (SDE) and conceptually different from an ordinary
(ODE) or a partial differential equation (PDE). σ is the amplitude of the noise
term and η describes Gaussian white noise which obeys:
1. hη(t)i = 0
2. hη(t)η(t′ )i = 2δ (t − t′ )
40
The formal solution is given by:
!
t
σ
ˆ
−t/t0 s/t0
v(t) = e v0 + ds e η(s)
0 m
as one can check easily by insertion into the Langevin equation. Here t0 = m/ξ
is the characteristic relaxation time of the system.
Obviously v is defined only through its averages, like the noise itself:
hv(t)i = v0 e−t/t0
′ 2 ′ t t′ s+s′
σ
ˆ ˆ
− t+t − t+t
′
v(t)v(t ) = v02 e t0 + e t 0 ds ds′ e t0 2δ s − s′
m
|0 0
{z }
t<t′ ´ t
= 0 ds 2e2s/t0 =t0 (e2t/t0 −1)
!
− t+t
′
σ2 σ 2 (t′ −t)/t0
=e t
0 vo2 − + e
mξ mξ
| {z }
=0 for t,t′ ≫t0
σ2 D E
⇒ v(t)2 =
mξ
Note that the linear terms in η have droped out and that the autocorrelation
decays exponentially, thus the system is well-behaved.
The equipartition theorem gives us:
1 D 2E 1
m v = kB T
2 2
⇒ ξv = ση(t) = ξ ẋ
1 t ′
ˆ
⇒ x(t) = x0 + dt ση(t′ )
ξ 0
⇒ hx(t)i = x0
1 t ′ t ′′ 2 ′
D E ˆ ˆ
(x(t) − x0 )2 = 2
dt dt 2σ δ t − t′′
ξ 0 0
1 2 !
= 2 2σ t = 2Dt
ξ
Here we identified the diffusion constant D from the one dimensional random
walk.
σ2 kB T
⇒ D= 2 = Einstein relation
ξ ξ
41
If we use for the friction coefficient Stoke’s law from hydrodynamics, ξ = 6πηR
with viscosity η we get:
kB T
⇒ D= Stokes-Einstein relation
6πηR
kB T 1/2
dx(t)
= −M ∇U + η (2.27)
dt ξ
U
= −D∇ + D1/2 η (2.28)
kB T
Here we have introduced the mobility M = 1/ξ. Due to the FDT, we only have
one relevant parameter, which we take to be D. The discretized version now
reads:
d
U
√
x(t + ∆t) = x(t) − D ∆t + 2D∆t N (0, 1) (2.29)
dx kB T
where N (0, 1) is the Gaussian distribution with vanishing mean and unit variance.
In order to be able to work with the standard Gaussian, we now explicitly use the
factor of 2 that above we have placed in the definition of the Gaussian white noise.
It is straight-forward to generalize this scheme to N particles with interaction
terms in U .
The BD-community did not yet converge on a few software packages and thus
there are many of them. Here are a few examples:
• HOOMD from Sharon Glotzer’s lab, also a MD-code with Langevin mode,
https://hoomd-blue.readthedocs.io
• Greens Function Reaction Dynamics (GFRD) from Pieter Rein ten Wolde,
event-based reactions based on exact solutions to the diffusion equation,
gfrd.org
42
• BrownDye from Gary Huber (Andrew McCammon group), uses APBS for
electrostatics, browndye.ucsd.edu
• MCell from Terry Sejnowski, Monte Carlo, often used in the neurosciences,
www.mcell.cnl.salk.edu
We finally note that in BD one deals with an effective solvent that in principle
should also mediate hydrodynamic interactions. The importance of hydrodynam-
ics for self-diffusion and molecular interactions is somehow debated and there are
many approaches to address this important issue. Because biomolecules and cells
are so small, they have a small Reynolds number:
ρvL
Re = ≪1 (2.30)
η
if we insert typical values for density ρ and viscosity η of water as well as for
velocity v and size L. Thus viscous forces dominate over inertial ones and one
has to solve the Stokes equation rather than the Navier-Stokes equation. The
most common approaches to hydrodynamics in biological systems are
For a recent review, compare Ulf Schiller et al, Mesoscopic modelling and simu-
lation of soft matter, Soft Matter 2018.
43
44
Chapter 3
We already discussed that polarizability and thermal rotation weakens the elec-
trostatic interaction between two molecules in the cell. We now consider this issue
further for a charged object (an ion, a biomolecule or an assembly of biomolecules)
immersed in a sea of other charges. These surrounding charges are counter-ions
(balancing the charge on the object) that can be complemented by co-ions (salt).
The fact that the counter- and co-ions are highly mobile and under permanent
thermal motion creates a cloud of charge distribution that screens and further
decreases the electrostatic interaction between the large objects.
45
Figure 3.1: Left panel: Polyelectrolyte such as a DNA molecule. Per base pair (bp,
distance between bp a = 3.4 Å), the DNA carries a charge of 2 e− and hence a line charge
2e−
density of λ = 3.4 Å
(image of DNA molecule taken from Wikipedia). Right panel:
Negatively charged head groups of the fatty acids in the plasma membrane resulting in
e−
an area charge density of σ = nm 2.
relevant biophysical questions are very different from the case of DNA. Because
of the barrier function of the lipid bilayer, we have to ask how charges arrange
themselves in its vicinity and how they can cross the bilayer. Obviously the dis-
tribution around charged lines and surfaces must be very different for geometrical
reasons.
As already mentioned in the beginning of this chapter, a DNA molecule can
2e−
be seen as a charged line with a linear charge density of λ = 3.4Å . To simplify
matters, we assume the DNA molecule to be an infinitely long straight line. Then
the DNA exhibits a cylindrical symmetry (compare figure 3.2a) and Gauss law
can easily be applied to determine the radial component of the electrostatic field:
λL
Gauss’ law: Er · 2πrL =
ǫ0 ǫ
λ
Electrostatic field: Er = (3.1)
2πǫ0 ǫr
´r
Φ = − a Er dr !
Electrostatic potential: λ r (3.2)
= − ln
2πǫ0 ǫ a
where a microscopic limit a was employed. The logarithmic electrostatic potential
Φ in equation 3.2 diverges for r → ∞. Thus, the boundary condition Φ(∞) = 0
cannot be used. One often encounters this logarithmic behavior in 2D systems.
For example, this means that one cannot calculate a simple formula for the flow
of a fluid around a cylinder (as one can for the flow of a fluid around a sphere in
3D).
For the straight line charge the same result as in equations 3.2 and 3.2 can be
obtained from the direct integration of Coulomb’s law (equation 1.56) or from
46
Figure 3.2: Cylindrical symmetry of a.) an infinitely long charged line and b.) an charged
plane with infinite surface.
Aσ
Gauss’ law: Ez · 2A =
ǫ0 ǫ
σ
Electrostatic field : Ez = (3.3)
2ǫ0 ǫ
σz
Electrostatic potential: Φ=− (3.4)
2ǫ0 ǫ
Two comments can be made concerning the results in equations 3.3 and 3.4.
Firstly, Φ increases linearly with the distance from the charged plane. Secondly,
the electric field jumps by σ/(ǫ0 ǫ) across the charged plane and does not de-
pend on the distance. As before, the same results can be obtained from explicit
integration or from solving the Poisson equation.
Besides its function as a diffusion barrier, the biomembrane can act as a parallel
plate capacitor (compare figure 3.3) if charges are separated to its both sides by
active processes such as ion pumps and transporters. Then we are actually dealing
with two oppositvely charged planes with electric fields according to equation 3.3:
σ
E+ = E− = (3.5)
2ǫ0 ǫ
47
Figure 3.3: Electrostatic potential (Φ, labeled in red) and ion concentration (c, labeled
in green) across the plasma membrane which is modeled as a parallel plate capacitor.
In the inter-membrane region, the potential decreases linearly whereas the concentration
follows ∼ e−q∆Φ/(kB T ) (equation 3.10). Since the electrical field of a charged plane does
not depend on the distance from the plane (compare equation 3.3), the net field outside
the membrane vanishes. Hence, there is a force on a charged test particle only if it is
within the lipid bilayer.
Outside the plasma membrane, E+ and E− cancel each other, whereas within
the membrane they add up:
σ
Electrostatic field: Einside = E+ + E− = (3.6)
ǫ0 ǫ
ˆ d
σ σd
Electrostatic potential difference: ∆Φ = − dz =− (3.7)
0 ǫ0 ǫ ǫ0ǫ
Q Aσ Aǫ0 ǫ
C= = = (3.8)
U |∆Φ| d
C ǫ0 ǫ µF
= ≈
A d cm2
where F denotes the physical unit Farad. This value membrane capacity has
µF
been measured experimentally. Moreover, the measure of 1 cm 2 is universally
used in experiments to determine the area of any given cell or membrane patch.
The concept of the biomembrane being a parallel circuit of a capacitor and an
ohmic resistance forms the basis of electrophysiology (theory of action potentials
according to Hodgkin and Huxley).
We now consider a single species of mobile ions that is free to distribute along z
(e.g. ions diffusing through the hydrophobic part or through an ion channel of the
membrane). At finite temperature T , there is a competition between electrostatic
forces and translational entropy. The concepts of energy stored in the form of
48
intracellular (mM) extracellular (mM) Nernst potential (mV)
K+ 155 4 -98
N a+ 12 145 67
Cl− 4 120 -90
Ca2+ 10−4 1.5 130
Table 3.1: Nernst potentials for some important ions in a typical mammalian muscle
cell. Because the Nernst potentials of the different ion species differ strongly, this ionic
distribution is an out-of-equilibrium situation. Resting potentials of excitable cells are
in the range of −50 to −90 mV .
Equation 3.10 was first formulated by the German physical chemist Walter Nernst
who won the Nobel prize in chemistry in 1920. It can be seen as Boltzmann’s
law for charges in an electrostatic potential (compare figure 3.3). In table 3.1
we give experimentally measured values for ion concentrations in a muscle cell.
The corresponding Nernst potentials are calculated in the last column. One sees
that they differ widely, proving that the distributions are out off equilibrium (ion
pumps and channels redistribute them against the thermal forces).
Our discussion showed that mobile charges will lead to concentration profiles
that depend on temperature and electrostatic potential. Therefore we now turn
to "electrostatistics", the field that combines these two elements.
1. high T or small charge density σ (no salt): In this case mean field theory
(MFT) can be used to derive the Poisson-Boltzmann theory.
49
Figure 3.4: Concentration profile of counter-ions (here: positive) and co-ions (here: neg-
ative) in a solution (a) without salt and (b) with salt, at a distance z from the charged
wall. For (b), the physiological concentration of additional salt is for instance in the
range of cs = 100 mM (e.g. NaCl).
50
Because we focus on the effect of a wall, we now rescale all distances with µ:
H X Ξ X
= + z̄i (3.13)
kB T r̄
i<j ij i
where
Z 2 lB Z 3 e4 n2d
Ξ= = 2πZ 3 lB
2n
2d = coupling strength (3.14)
µ 8π(ǫ0 ǫkB T )2
In equation 3.13, we rescaled the system such that only one dimensionless pa-
rameter, namely the coupling strength (equation 3.14), determines the behavior
of the whole system.
At this point, we managed to end up with only one dimensionless parameter,
that defines two asymptotic limits of interest2 :
1. Ξ ≪ 1: This is the case if the system has a low charge density, a low
valency and/or a high temperature. One can perform an expansion in small
Ξ (mean-field theory) and ends up with the Poisson-Boltzmann theory.
In order to understand the difference better between the two limits, we use charge
neutrality
Ze
= σ = en2d (3.15)
πa2⊥
This shows that Ξ determines the ratio between a⊥ and µ (compare figure 3.5).
For Ξ ≪ 1, the lateral distance between the counter-ions is smaller than their
average distance from the wall and they form a 3D cloud that has no structure in
the lateral direction; therefore a mean field theory in z-direction is sufficient. For
Ξ ≫ 1, the lateral distance between the counter-ions is larger than their average
distance from the wall and they form a 2D layer on the wall. For very strong
coupling, this condensate can become a crystal.
2
A complete treatment can be given using statistical field theory. For more details, compare
the overview in R.R. Netz, European Physical Journal E 5: 557, 2001.
51
Figure 3.5: The two complementary limits considered here. (a) In the high temperature
limit, the lateral distance of the counter-ions is smaller than the vertical extension and
thus we get a 3D fluid. This situation is described by a mean field theory in z-direction
(Poisson-Boltzmann theory). (b) In the low temperature limit, the lateral distance is
larger than the vertical extension and we get a 2D condensate (strong coupling limit).
This results in
e
⇒ ∆Φ = − · n0 · exp − keΦ
BT
Poisson-Boltzmann equation
ǫ0 ǫ
(3.18)
The Poisson-Boltzmann equation (PBE) is a non-linear differential equation of
second order which is in general hard to solve analytically. In MD simulations,
one usually employs PB-solvers (e.g. DelPhi, APBS, MIBPB, etc). There are
only few cases for which it can be solved analytically.
Luckily, this is the case for the example of the charged wall. The boundary
conditions are given by the charge neutrality of the whole system and by |E(∞)| =
| − Φ′ (∞)| = 0:
´∞ ´∞
σ = − 0 dz ρ(z) = ǫ0 ǫ 0 Φ′′ dz
|{z}
charge density
of counter-ions
= ǫ0 ǫ( Φ′ |z=∞ − Φ′ |z=0 )
| {z }
= −E(∞) = 0
σ
⇒ Φ′ |z=0 = −
ǫ0 ǫ
52
point charge charged wall with counter-ions (PBT) with salt (DH)
Φ 1/r z ln(z) exp(−κz)
E 1/r2 const 1/z exp(−κz)
Table 3.2: Distance dependence of the electrostatic potential Φ and the electrostatic field
E for different systems. Note that in comparison to a point charge a spatially extended
distribution like the charged wall strengthens the interaction, whereas the presence of
counter-ions (Poisson-Boltzmann theory) weakens the interactions. If, in addition, salt
is added to the solution, the interaction is weakened to an even higher extent.
With the boundary conditions, we get the analytical solution for the charged
wall:
2kB T
z+µ
Electrostatic potential Φ(z) = ln + Φ0 (3.19)
e µ
1 1
Counter-ion density n(z) = · (3.20)
2πlB (z + µ)2
Recall, that without counter-ions Φ ∼ z and E = const (compare equation 3.2
and equation 3.3; compare also table 3.2). This is now changed to a logarithmic
scaling of the potential since a cloud of counter-ions surrounds any charged object
and thus weakens the electrostatic potential. In other words, the charged wall is
"screened" by the counter-ions. Together with the cloud or layer of counter-ions,
the charged wall forms an electrostatic "double layer".
53
lB µ lDH
7Å 1 nm 1 nm
Table 3.3: Values for the three electrostatic length scales Bjerrum length lB , Gouy-
Chapman length µ and Debye-Hückel screening length lDH at physiological conditions.
Note that the three electrostatic lengths are very similar (all around 1 nm).
54
QV
1. No salt added to the solution, hence κ → 0 ⇒ Φ = . This limit
4πǫ0 ǫr
results in the well-known Coulomb law.
2. Point charge (R → 0): Then the potential takes the form
QV · exp(−κr)
Φ= Yukawa potential (3.28)
4πǫ0 ǫr
Equation 3.28 is the Green’s function (or propagator) for the linear Debye-
Hückel theory. Like for the Coulomb interaction, one can calculate Φ for
any extended object (i.e. a line, a plane, etc.) by superposition of the
propagator.
55
Figure 3.6: The formation of river deltas as a consequence of the interaction of two
charged particles in salty solution. (a) A river flows from the mountains into the sea. On
its way, negatively charged silica particles dissolve in the water which repel each other
due to their charge. As the low-salt water of the river meets the sea water with high
salinity, the repulsion between the particles is screened. They aggregate and, hence, form
the river delta. (b) Effect of the salt ions on the potential energy. Screening lowers the
energy barrier responsible for the repulsion. The corresponding description is known as
DLVO-theory in colloidal sciences.
Figure 3.7: (a) PB solution for the potential Φ ∼ ln cos2 (K · z) (red) and the counter-
ion density n ∼ 1/ cos2 (K · z) (blue) between two charged walls with charge densities
σ1 = σ2 = σ. (b) The human knee is stabilized by cartilage containing hyaluronic acid
(HA). Hyaluronic acid is a long, high molecular mass polymer of disaccharids, which
is negatively charged and therefore responsible for the disjoining pressure caused by its
counter-ions.
56
The two boundary conditions are:
symmetry : Φ′ (0) = 0
ˆ d/2 d/2
d
ˆ
′′ ′
charge neutrality: σ= ρ dz = ǫ0 ǫ Φ dz = ǫ0 ǫΦ
0 0 2
d σ
⇒ Φ′ =
2 ǫ0 ǫ
This results in the exact solution of the PBE:
kB T h 2 i
Potential Φ(z) = ln cos (K · z) (3.30)
e
n0
Counter-ion density n(z) = 2
(3.31)
cos (K · z)
where n0 denotes the counter-ion density at the mid plane and k denotes a con-
stant which follows from the boundary condition:
!
2kB T · K
′ d σ K ·d
Φ = = tan (3.32)
2 ǫ0 ǫ e 2
p = kB T · n0 = 17 atm (3.33)
where 1 atm ∼ = 105 P a. It can be seen directly that the disjoining pressure is
very large and this has many applications in biological systems. For example,
disjoining pressure can be found in joints and is actually the reason why we can
go jogging (compare figure 3.7b).
57
Figure 3.8: (a) Phase diagram showing regions of attraction and repulsion as a function
of plate separation d/µ and coupling strength Ξ. Original phase diagram by Moreira &
Netz: "Binding of similarly charged plates with counter-ions only." Phys Rev Lett, 87(7):
078301, Aug 2001. (b) Small strip of two charged walls with an in-between counter-ion.
simply a superposition of the solution of two single charged walls (equation 3.25).
One then gets
Thus, the DH solution for Φ (as well as for n and p, respectively) decays expo-
nentially with the distance and, hence, the interaction is short-ranged.
The counter-ion density between two charged walls in the strong coupling limit
turns out to be relatively flat. In detail it is constant in cero order and parabolic
in first order of the virial expansion. Thus superficially it appears to be similar
to the PB result. In practise, however, the results are very different, because
one finds that the two equally charged walls can in principle attract each other.
Whether the interaction between the two planes is attractive or repulsive depends
on the distance d and the coupling strength Ξ, as shown in the phase diagram in
figure 3.8a. The very fact that attraction can occur offers a solution to our DNA
riddle.
A simple explanation for this behavior can be given as follows: consider the
condensed situation as sketched in figure 3.8b. Because the counter-ions condense
with a relatively large lateral distance to each other, we neglect their interaction
and only consider the interactions of one counter-ion with the wall in a small strip
q
with area A = 2σ . There are three contributions to the electrostatic energy now:
the two interactions of the counter-ion with the two walls and the interaction of
the walls with each other:
Uel
= 2π(lB /e2 )qσx + 2π(lB /e2 )qσ(d − x) − 2π(lB /e2 )σ(σ · A)d
kB T
= π(lB /e2 )σqd = 2π(lB /e2 )σ 2 Ad (3.35)
The energy is minimal for d → 0 which leads to attraction of the two charged
58
walls. For the electrostatic and the entropic pressure we get
∂ Uel −2πlB σ 2 kB T
electrostatic pressure: pel = − =
∂d A e2
kB T 2σkB T
entropic pressure: pen = =
A·d qd
e2
⇒ balanced at equilibrium distance d= (3.36)
πlB qσ
The strong coupling limit is biologically relevant, because for n2d = 1 nm−2 it
can be reached with trivalent counter-ions. In fact, the charged polymer DNA
uses many multivalent counter-ions such as speridine and spermine which support
DNA condensation in the nucleus. Again the existence of an equilibrium distance
also has consequences in other sciences. E.g. it explains why clay particles can
be swollen only to a certain distance.
59
Figure 3.9: Chromatin has a highly complex structure with several levels of organization.
Source: Pierce, Benjamin. Genetics: A Conceptual Approach, 2nd ed. (New York: W.
H. Freeman and Company 2005).
DNA− Henderson-Hasselbalch
⇒ pK = pH − log10 (3.40)
[DNA] equation
The
pK corresponds to the pH at which half of the groups have dissociated
( DNA− = [DNA]).
For DNA, we find pK = 1 which implies that DNA is a very strong acid. In cells,
pH = 7.34. With the Henderson-Hasselbalch equation the fraction of dissociated
DNA can immediately be calculated.
DNA−
= 106.34
[DNA]
Thus, DNA in the cell is completely dissociated and, therefore, carries a line
charge density of
Now we want to focus on DNA packing in viruses. Actually, a virus is not a living
object per definition, but rather genetic material, i.e. DNA or RNA, packed into
a protein shell, the so-called capsid. Typically, the diameter of a capsid is in the
range of 10s to 100s of nm. Some viruses, e.g. HIV, are in addition wrapped by
a lipid bilayer (and are then called "enveloped virus").
60
Figure 3.10: Simple model of viral DNA packed in the capsid of the φ29 bacteriophage.
As we shall see in the following, the RNA and DNA, respectively, in viruses is
very densely packed. Take for instance the φ29 bacteriophage (a virus infecting
E.Coli): Its capsid can be approximated as a sphere of radius Rcapsid = 20 nm
containing 20 kbp (corresponding to L = 2 · 104 · 0.34 nm = 7 µm) DNA. We
assume Vbp ≈ 1 nm3 (compare figure 3.10). The packing ratio in the capsid can
be computed directly:
2 · 104 nm3
≈ 0.6 (3.42)
4π
3 (20 nm)3
Comparing this value with the maximal packing density of spherical objects into a
crystal (≈ 0.71) it can be concluded that DNA packed into a viral capsid must be
close to a crystalline structure. Indeed this can be shown by electron microscopy.
If we now pack DNA with the line charge density λ = 2e/(3.4 Å) into the virus,
how much electrostatic energy do we have to put into the system? Electrostatic
energy is the work to bring a charge distribution into its own field and is known
to be
1
Uel = Φ(~r) · ρ(~r) dr (3.43)
´
2
1 1 4π 3
Φ(r) = · · ρ· r (3.44)
4πǫ0 ǫ r | {z3 }
| {z }
point charge in origin charge in smaller sphere
61
Figure 3.11: Left panel: Experimental set-up of the portal motor force experiment. A
single φ29 packaging complex is tethered between two microspheres. Optical tweezers
are used to trap one microsphere and measure the forces acting on it, while the other
microsphere is held by a micropipette. Right panel: Internal force as the function of %
genome packed. Note that 100% of the genome corresponds to 7 µm DNA and that the
work is obtained by integration of the force (grey area). Images and caption text (partly)
taken from reference in footnote 5.
For the total work, we have to add up shell after shell of the sphere:
ˆ R
⇒ Uel = dr Φ(r) · ρ · 4πr2
0
R 4π 2 4 4π 2 5
ˆ
= dr ρ r = ρ R
0 3ǫ0 ǫ 15ǫ0 ǫ
1 3Q2
= · (3.45)
4πǫ0 ǫ 5R
Here the factor 1/2 does not arise because every contact is counted only once as
we gradually build up the sphere. For our example of the φ29 bacteriophage, we
have Q = 2e/bp · 20kbp and hence
The work needed to pack the DNA into the viral capsid has been measured
in a single molecule experiment5 as sketched in figure 3.11. However, in this
experiment the work was determined to be much smaller than the one estimated
above:
1
Wexp ≈ 7000 nm 60 pN = 2.1 · 105 pN · nm (3.47)
2
Obviously the above estimate was much too high because we neglected the effect
of the counter-ions.
5
DE Smith et al. "The bacteriophage straight φ29 portal motor can package DNA against a
large internal force." Nature, 413(6857):748–52, Oct 2001.
62
There are N = 4 · 104 counter-ions packed with the genome (corresponding to
2 counter-ions/bp). We now assume complete neutralization of the charges and
consider only the loss of entropy due to the DNA volume:
Vf ree
Uci = N kB T · ln (3.48)
Vcapsid
The volume Vf ree is that of the screening cloud (recall that L ≈ 7 µm, RDN A ≈
1 nm, Rcapsid ≈ 20 nm and lDH ≈ 1 nm).
h i
Vf ree = Lπ (RDN A + lDH )2 − RDN
2 4
A = 6.6 · 10 nm
3
(3.49)
4π 3 2 4 3
Vcapsid = R − LπRDN A = 1.2 · 10 nm (3.50)
3 capsid
⇒ Uci = 3 · 105 pN · nm (3.51)
This result is much closer to the experimental value. A full analysis had to also
include the effect of bending the DNA, which requires polymer physics.
Finally the pressure inside the capsid can be calculated:
N kB T 4 · 104 · 4.1 pN · nm
p = = 4π 3
3 (20 nm)
V
pN
≈ 5 = 50 atm (3.52)
nm2
This is a huge counter-ion pressure inside the capsid, as was also experimentally
confirmed.
63
(a) (b)
(c) (d)
Figure 3.12: (a) Molecular rendering of the structure of the capsid of hepatitis B virus
(HBV). (b) Assembly isotherms at different salt concentrations. (c) Fit of equilibrium
constant as a function of temperature and salt concentration. (d) Scaling of equilibrium
constant with salt concentration. Experimental data from Slotnick group, theory by van
der Schoot group.
For HBV, K has been measured as a function of temperature T and salt con-
centration cs 6 , compare figure 3.12b-d. Because this virus assembly from 120
capsomers in an all-or-nothing manner, we do not have to consider intermediates
and can write a law of mass action:
[capsid]
K= (3.55)
[capsomer]120
64
In a theoretical analysis, it has been shown that these experimental results can
be fitted nicely using Debye-Hückel theory 7 . We start from the surface potential
of a sphere of radius R in Debye-Hückel theory (equation 3.27):
QV QV lDH
φR = = . (3.56)
4πǫ0 ǫ · (1 + κR)R 4πǫ0 ǫ R(R + lDH )
With lDH = 1nm and R = 14nm we can write (R + lDH ) ≈ R. Therefore our
final result for the salt-dependent part of the equilibrium constant reads
2
∆Gelectro 1 Q lDH lB
ln K = − = − kB T (3.58)
kB T 2 e R2
Thus ln K should scale linearly with the screening length lDH and therefore with
−1/2
cs , exactly as it is observed experimentally, compare figure 3.12d.
7
W.K. Kegel and P. van der Schoot, Biophysical Journal 86: 3905, 2004
65
66
Chapter 4
In a cell, lipid bilayers partition space into functional compartments. This central
aspect of lipid bilayers must have been crucial for the development of life. Lipid
bilayers are the carriers of many vital processes, including ion separation and
transport as well as protein activity. In general, cell membranes regulate the
transfer of material and information in and out of cells.
Due to its low bending energy and the thermal environment, the lipid bilayer
is in continuous motion. In order to describe the energetics and statistics of
membranes, we have to introduce a mathematical description of surfaces and then
to identify the corresponding energy (Helfrich bending Hamiltonian). Therefore,
we start with a crash course in differential geometry1 . We then discuss the Helfrih
bending Hamiltonian in much detail and its consequences for shapes of minimal
energy and thermal fluctuations around these shapes. As a reference point, we
always discuss surfaces under tension (e.g. soap bubbles or oil droplets). Finally
we discuss the physics of red blood cells, whose shapes and fluctuations can be
described well by surface Hamiltonians. However, in contrast to pure membranes,
the presence of the actin-spectrin network makes it necessary to add additional
terms to the interface Hamiltonian.
1
There are many books on differential geometry, for example the one by Michael Spivak
(Comprehensive introduction to differential geometry, vols 1-5, 1979). Here are two books in
German that are especially helpful for membrane physics: MP do Carmo, Differentialgeometrie
von Kurven und Flächen, 3rd edition Vieweg 1993; JH Eschenburg and J Jost, Differentialge-
ometrie und Minimalflächen, 2nd edition Springer 2007. The classical review on vesicle shapes
is Udo Seifert, Configurations of fluid membranes and vesicles, Advances in Physics 46: 13-137,
1997. A great resource is also the script by JL van Hemmen, Theoretische Membranphysik:
vom Formenreichtum der Vesikel, TU Munich 2001, available at http://www.t35.physik.tu-
muenchen.de/addons/publications/Hemmen-2001.pdf from the internet.
67
Figure 4.1: Membranes and polymers can be mathematically described as two-
dimensional surfaces and one-dimensional curves, respectively, in a three-dimensional
space.
Consider a curve in 3 dimensions, e.g. a helical curve with radius R and pitch
z0 = b · 2π
ω (figure 4.2a), parametrized by an internal coordinate t:
x1 (t) R · cos(ωt)
~r(t) = x2 (t) = R · sin(ωt) (4.1)
x3 (t) b·t
In the limit b → 0, the helix becomes a circle, and in the limit b → ∞, it becomes
a straight line. For the velocity of the helix we get:
d~r −Rω · sin(ωt)
~v = ˙
= ~r = Rω · cos(ωt)
dt
b
p
⇒v = R2 ω 2 + b2 > Rω
An important quantity when describing a curve is its arc length L which is inde-
pendent of the parametrization that was chosen.
t1
t=t(u) ˆ t1 d~
r du
ˆ
˙
L = dt ~r = dt ·
t0 t0 du dt
ˆ u1
d~
r
= du (4.2)
u0 du
68
Figure 4.2: a.) Helical curve with radius R and pitch zo = b·2π/ω. The tangential vector
~t, the normal vector ~n and the binormal vector ~b are sketched in blue. b.) "Kissing circle"
at a point P (s) with radius R(s) = κ−1 (s).
ˆ t
s(t) = dt′ ~r˙ (t′ ) (4.3)
t0
can be used to parametrize the curve since it increases strictly with t (ṡ = ~r˙ > 0)
parametrization by
⇒ r = r(s) = r(t(s)) (4.4)
arc length (PARC)
√ ds
v= R 2 ω 2 + b2 =
= const
dt
R · cos( ωs
v )
s
⇒ s = v · t ⇒ t = ⇒ ~r = R · sin( wsv )
v bs
v
69
The co-moving frame
The co-moving frame (also called "Frenet frame") of a curve consists of three
mutually perpendicular unit vectors:
~r˙ PARC
d~
r
· ds
~t(s) := ds dt d~r
tangential vector = = (4.5)
|~r˙ | ds
dt ds
−1
d~t d~t 1 d~t
normal vector ~n(s) := · = (4.6)
ds ds κ ds
binormal vector ~b(s) := ~t(s) × ~n(s) (4.7)
The normalization in equation 4.5 is not required for PARC. Therefore, PARC is
also called the "natural parametrization". In equation 4.6 we defined the curva-
ture κ of the curve at a given point
d~t
:= (4.8)
κ
ds
which defines the radius of curvature R(s) = κ−1 . This is the radius of the
so-called "kissing circle" at that specific point (figure 4.2b).
E.g. for the helical path we get:
− Rω · sin( ωs
v )
d~r(s) Rωv ωs R 2 ω 2 b2 √
~t = = v · cos( v ) ⇒ |~t| =
+ 2=1
ds v 2 v
b
v
Rω 2 Rω 2
− Rω
2
· cos( ωs κ = =
d~t v22 v ) v2 R2 ω 2 +b2
= − Rω · sin( ws
⇒ 1
ds v2 v ) = 1 <
0 2
R+ b 2
Rω
R
The curvature of the helical path is smaller than for a circle. In the limit b → 0,
we have κ = 1/R, denoting a perfect circle. In the limit b → ∞, κ vanishes,
denoting a straight line.
The derivatives of the vectors of the co-moving frame are described in the same
basis through the Frenet formulae:
d~t
= κ~n
ds
d~n
= −κ~t +τ~b (4.9)
ds
d~b
= −τ~n
ds
where we introduced the torsion τ :
d~b d~n ~
τ = − · ~n = ·b (4.10)
ds ds
70
Figure 4.3: a.) Surface in a 3D space with tangential vectors ∂x f~ and ∂y f~ and unit
normal vector ~n perpendicular to the tangential vectors. b.) Plane (yellow) containing ~n
and rotating. c.) For each position Θ of the rotating plane a curvature can be determined.
τ measures how strongly the curve is twisted out of the plane. E.g. for the helical
path
d~b bω bω b→∞
τ = −~n · = 2= 2 2 −−−−→ 0
ds v R ω + b2 or b→0
4.1.2 Surfaces in 3D
The tangential vectors ∂x f~ and ∂y f~ span the tangential plane (compare figure
4.3a). The unit normal vector then is defined as
∂x f~ × ∂y f~
~n =
~
(4.12)
∂x f × ∂y f~
Note that in contrast to the case of space curves, we do not normalize the tan-
gential vectors.
In order to introduce definitions for the curvature, we can construct a plane
containing ~n which we then rotate by 180 degrees through a given point (x, y)
on the surface, as sketched in figure 4.3b. The kissing circle for each span curve
defined by a rotation angle Θ of the plane gives us a curvature in this certain
direction (figure 4.3c).
The curvature will have a minumum κ1 and a maximum κ2 , the so-called "principal
curvatures". With these two curvatures and the radii of the corresponding kiss-
71
ing circles R1 and R2 , respectively, we can define two important concepts:
!
κ 1 + κ2 1 1 1
Mean curvature: H := = + (4.13)
2 2 R1 R2
1
Gaussian curvature: K := κ1 · κ2 = (4.14)
R1 · R2
Example
Radii of
R1 = R, R2 = ∞
kissing R1 = R2 = R R1 = −R2
(straight line)
circles
Mean
1 1
curva- H= R H=0 H= 2R
ture
Gaussian K = R12 > 0 (cannot
K = − R12 < 0 K = 0 (can be
curva- be mapped onto
1 mapped onto plane)
ture plane)
Table 4.1: H and K can be used to classify surfaces. For the examples shown here,
K is constant, and hence each point on the surface is elliptic, hyperbolic or parabolic,
respectively.
72
where χ is the so-called "Euler characteristic". It can be used to calculate
the number of handles G ("genus") of a surface:
χ = 2−2·G (4.16)
fA (x, y) dxdy and K = K(x, y). To this end, we calculate the three 2 × 2 ma-
trices g, h and a. Let us first define the symmetric matrix g, also called "first
fundamental form" or "metric tensor":
~2
∂x f~ · ∂y f~
∂x f
gij := ∂i f~ · ∂j f~ = metric tensor (4.18)
~2
∂x f~ · ∂y f~ ∂y f
2
1 ∂y f~ −∂x f~ · ∂y f~
−1
⇒ gij = ~2
(4.19)
det g −∂x f~ · ∂y f~ ∂x f
2
where det g = ∂x f~ × ∂y f~ . gij depends on ∂x f~ and ∂y f~, but not on the unit
which, in contrast to the metric tensor, depends on the unit normal ~n.
With the matrices g and h, the matrix a can be defined:
Like curvature κ and torsion τ tell us how the normal changes along a space
curve, the Weingarten matrix tells us how the normal changes along a surface:
aij ∂j f~
P
∂i~n = − j (4.23)
73
Euler characteristic Genus
Topological 2−χ
Object χ = F −E+V G =
equivalence 2
Cube
χ = 6 − 12 + 8
G = 0
a.) Sphere = 2
Tetrahedron
χ = 4−6+4
G = 0
= 2
punctured cube
χ = 16 − 32 + 16
c.) Torus G = 1
= 0
d.) Double
χ = −2 G = 2
torus
Table 4.2: a.) The sphere is topologically equivalent to a cube and a tetrahedron,
respectively. χ can also be calculated from the Gauss-Bonnet theorem (equation 4.15):
dA K = 4πR2 · R12 = 2π·2. b.) The Euler characteristic is additive over multiple bodies.
¸
The more bodies, the larger χ and the more negative G. c.) The torus is topologically
equivalent to toroidal polyhedra, e.g. a punctured cube. Note that G = 1, denoting
that the object has one handle. d.) For topologically more complex structures, like the
double or triple torus, it is more reasonable to determine the number of handles G and
then calculate χ from equation 4.16 than to find a topologically equivalent polyhedron.
Generally we find: The more handles, the larger G and the more negative χ.
74
(a) Plane (b) Cylinder
Figure 4.4
From the Weingarten matrix we now can compute the mean curvature H and
the Gaussian curvature K for any given surface f~ (without proof):
det h
K = det a = det g (4.24)
H = 12 tr a (4.25)
In the following, we will use this recipe for some important examples.
x 1 0 0
~ ~ ~
f = y ⇒ ∂x f = 0 , ∂y f = 1 , ~n = 0
0 0 0 1
!
1 0 p
⇒g= ⇒ det g = 1, dA = det g dxdy = dxdy
0 1
!
0 0
h=a= ⇒ H=K=0
0 0
75
2. Cylinder (figure 4.4b, internal coordinates φ and z)
R · cos φ −R · sin φ 0 cos φ
f~ = R · sin φ ⇒ ∂φ f~ = R · cos φ , ∂z f~ = 0 , ~n = sin φ
z 0 1 0
!
R2 0
⇒g= ⇒ det g = R2 , dA = R dφdz
0 1
! !
−R 0 − R1 0
h= a=
0 0 0 0
1
⇒ H=− , K=0
2R
x 1 0
~
f =
y
~ ~
⇒ ∂x f = 0 , ∂y f = 1 ,
h(x, y) ∂x h ∂y h
1 ∂x h
~n = q
∂y h
2 2
1 + (∂x h) + (∂y h) 1
!
1 + (∂x h)2 ∂x h · ∂y h
⇒g= ⇒ det g ≈ 1 + (∂x h)2 + (∂y h)2 = 1 + (∇h)
~ 2
∂x h · ∂y h 1 + (∂y h)2
76
Figure 4.5: Sketch of a bent lipid double layer.
1 ~ 2
q
dA ≈ ~ 2 dxdy ≈ [
1 + (∇h) 1 + (∇h) ] dxdy (4.26)
| 2 {z }
|{z}
reference excess area, >0
plane
!
1 ∂xx h ∂xy h
h= √ ≈a
det g ∂xy h ∂yy h
1 1
H ≈ (∂xx h + ∂yy h) = ∆h (4.27)
2 2
K ≈ ∂xx h · ∂yy h − (∂xy h)2 (4.28)
ˆ
H = dA f (H, K)
77
One can expand the Hamiltonian in small curvatures (or, in other words, in small
1/R) up to order 1/R2 to obtain2
Helfrich-Canham
H= dA {σ + 2κ(H − c0 )2 + κ̄K
´
Hamiltonian
Recall that H ∼ O(1/R) and K ∼ O(1/R2 ). Higher order terms have also been
investigated but lead to very complex structures.
The Helfrich-Canham Hamiltonian contains four material parameters:
78
(a)
(b)
Figure 4.6: Spontaneous curvature of the membrane due to asymmetries caused by a.)
embedded proteins (adsorbed proteins would have similar effect) and by b.) different
lipid compositions in outer versus inner leaflet.
This is why stability requires κ+ and κ− both to be larger than zero. This implies
−2κ < κ < 0, so κ is expected to be small and negative.
One can use the elasticity theory for linear isotropic material to derive the two
elastic moduli of the thin membrane as a function of the two elastic moduli of
the bulk material:
Ed3 −Ed3
κ= , κ = (4.29)
12(1 − ν)2 6(1 + ν)
where E is the Young’s modulus of the material and ν its Poisson’s ratio. The
bulk material relevant here is the hydrocarbon fraction of the bilayer. d = 4nm
is its thickness. With the value ν = 0.5 for incompressible material, a bending
79
Figure 4.7: The force to pull a tether out of a vesicle scales with the square root of the
tension. P. Basserau et al. Advances in Colloid and Interface Science 208 (2014) 47-57.
• Energetics deals with the question what is the surface with minimal en-
ergy. These surfaces have to solve the Euler-Lagrange equations, also called
shape equations δH δ
δh = 0. Here δh is a functional derivative.
which is a path integral or functional integral (integral over all possible func-
tions).
Together the minimal energy state and the fluctuations around it describe the
main features of interest.
80
where L is the length of the cylinder and F the pulling force. We minimize for
R to get: r
κ
R= (4.31)
2σ
At equilibrium, membrane energy and pulling energy should balance and thus
√
F = 2π 2κσ (4.32)
Therefore the force F scales as the square root of σ, as has been shown experi-
mentally, compare figure 4.7. In reverse, the force can be used to measure σ.
81
Figure 4.8: Parametrisation. The blue curve represents the part of the membrane that
adheres to the particle or is covered by the coat; the red line is the free part of the
membrane. Figure taken from Foret, Lionel. "Shape and energy of a membrane bud
induced by protein coats or viral protein assembly." The European Physical Journal E
37.5 (2014): 42.
Importantly, we can express r(s) and z(s) by means of the the tangential angle
ψ(s) as
ṙ = cos ψ(s) ż = − sin ψ(s) (4.35)
Our aim is to calculate the principal curvatures C1 and C2 . Hence, we calcu-
late the metric g, the normal vector ~n, the second fundamental form h and the
Weingarten matrix a. We find
!
1 0
g= (4.36)
0 r2
cos φ sin ψ
~n = sin φ sin ψ (4.37)
cos ψ
!
−ψ̇ 0
h= (4.38)
0 −r sin ψ
!
−ψ̇ 0
a= (4.39)
0 − sin ψ/r
Thus, the principal curvatures are given by the eigenvalues of a, C1 = −ψ̇ and
C2 = − sin ψ/r. Then the Helfrich Hamiltonian reads with the mean curvature
2H = C1 + C2 ( )
ψ 2
κ
ˆ
H = dsdφ σ + ψ̇ + sin r. (4.40)
2 r
82
because of the spherical geometry we do not need a second Lagrange multiplier,
as variations of the contour endpoints are independent. We define an action
H
ˆ ˆ
S[r(s), ψ(s)] = + dsγ(s)(ṙ − cos ψ) = dsL(ψ, ψ̇, r, ṙ) , (4.41)
2πκ
with a Lagrange function L
2
1 sin ψ r
L(ψ, ψ̇, r, ṙ) = ψ̇ + r+ + γ(ṙ − cos ψ) , (4.42)
2 r λ2
p
where λ = κ/σ defines the characteristic lengthscale of the membrane. The
Euler-Lagrange equations are the solutions to the variational problem
d ∂L ∂L
δS = 0 ↔ − , (4.43)
dt ∂ q̇k ∂qk
where qk = {r, ψ}. Thus,
83
4
α=0.15π α=0.28π α=0.42π α=0.57π α=0.71π α=0.85π
3 λ/R=0.1
2
0
-2 -1 0 1 2
4
3
λ/R=0.16
2
0
-2 -1 0 1 2
4
3
λ/R=0.32
2
0
-2 -1 0 1 2
4
3
λ/R=1
2
0
-2 -1 0 1 2
4
3
λ/R=3.2
2
0
-2 -1 0 1 2
6
4 λ/R=3.2
2
0
-6 -4 -2 0 2 4 6
6
4 λ/R=10
2
0
-6 -4 -2 0 2 4 6
Figure 4.9: Membrane shapes for different opening angles α and different values of λ/R.
In blue: the membrane that adheres to the particle or which is covered by the coat. In
red: the free part of the membrane. Figure taken from Foret, Lionel. "Shape and energy
of a membrane bud induced by protein coats or viral protein assembly." The European
Physical Journal E 37.5 (2014): 42.
We now aim for a phase diagram of particle uptake3 . For simplicity, now we
neglect the contribution from the free membrane. In general, for a particle to
be taken up, it has to be adhesive and the adhesion energy has to balance the
bending energy. The Helfrich Hamiltonian reads
ˆ
H = dA(2κH 2 + σ) − wAad (4.48)
3
The two classical papers on this subject are: R Lipowsky and HG Doebereiner, Vesicles in
contact with nanoparticles and colloids, Europhysics Letters 43:219-225, 1998; and M Deserno,
Elastic deformation of a fluid membrane upon colloidal binding, Phys Rev E 69:031903, 2004.
84
A B
Figure 4.10: (A) Particle uptake by membrane wrapping. (B) Phase diagram of particle
wrapping by membranes as function of surface tension σ and adhesion energy w. Taken
from M Deserno, PRE 2004.
where Aad is the adherent membrane area and w the adhesion energy area density.
A typical value would be w = 0.1 mJ/m2 .
We next consider a sphere of radius R. If we denote the angle α to describe where
the contact line between membrane and sphere is located (compare figure 4.10A),
then the wrapping variable z = 1 − cos α will run from 0 to 2 as the membrane
wraps the particle. If we neglect the contributions from the bending of the free
membrane, we get
E = 4πzκ + πR2 z 2 σ − 2πR2 zw (4.49)
The first term is the bending energy, which is independent of radius. The second
and last term are the surface tension and adhesion energy terms, respectively.
While both have the same R2 -scaling, they have different scaling with z. The
last term has the trivial z-scaling. The second term however scales as z 2 , because
here the excess area Aexcess = Aad − Aprojected matters:
E
Ē = = 4z + σ̄z 2 − w̄z = −(w̄ − 4)z + σ̄z 2 (4.51)
πκ
where σ̄ = σR2 /κ and w̄ = 2wR2 /κ. This energy function gives rise to a phase
diagram as shown in figure 4.10B. For w̄ < 4, the energy is always positive and
no wrapping can occur. This corresponds p to the free state. Note that w̄ = 4
translates into a minimal radius R = 2κ/w ≈ 50 nm, below which uptake is
not possible. For w̄ > 4 + 4σ̄, the minimal energy is found for z = 2, the fully
wrapped state. Inbetween there is a parameter region where the minimum lies
at a finite value of z, here the partially wrapped state is stable.
85
Figure 4.11: Schematics of a lipid vesicle with constant surface A and volume V . Al-
though lipid bilayers are permeable to water, there are always some ions caged in the
vesicle, fixing an osmotic pressure, which keeps the volume constant. Also the number of
lipids on the outside Nout and inside Nin is constant, because of the energy barrier that
does not allow for the lipids to jump from one side of the membrane to the other or to
escape.
In this section we are looking at closed surfaces; therefore κ is irrelevant due to the
Gauss-Bonnet theorem. Also c0 = 0 because we assume symmetric membranes.
We add a term −pV to control the volume. In practice, one prepares a suspension
of vesicles, e.g. by ultrasound or electoporation acting on a lipid-water mixture,
and then selects vesciles of interest, e.g. with optical tweezers or a micropipette.
Each vesicle then has fixed values for area and volume which can be measured
with e.g. video or fluorescence microscopy. Using A = 4πR02 , one can define the
radius of the equivalent sphere. Then the only relevant parameter of the system
is the reduced volume v:
V
v = 4π 3
3 R0
Each vesicle class has an individual value for v, and v ≤ 1 should be always
fulfilled; v < 1 describes a deflated vesicle with excess area for non-trivial shapes.
Shape with v = 1 is a sphere and has the optimal VA ratio. Note that
1
ˆ
H = 2κ dA H 2 = 2κ 4πR2 · = 8πκ = const.
|{z} R2
for sphere
which indicates that the solutions are independent of rescaling (this is part of a
more general property called conformal invariance).
In order to obtain a phase diagram as a function of v, we have to solve the
corresponding Euler-Lagrange equations (shape equations). These are derived by
varying the surface in normal direction
86
Figure 4.12: Catenoid as an example of a minimal surface, which is not compact and has
H = 0.
and then asking for the first ǫ-derivative of the energy functional to vanish for
arbitrary φ(x, y) (one can show that a tangential variation does not matter in
this order). The result can be written as4 :
p + 2σH − 2κ(2H H 2 − K − ∆H) = 0 Euler-Lagrange equation
(4.53)
where ∆ is the Laplace-Beltrami operator (only for the almost flat membrane
we get the Laplace operator ∆ = ∂ 2x + ∂ 2y ). The Euler-Lagrange equation is a
partial differential equation (PDE) of the 4th order with a famous limit for κ = 0,
namely the Laplace law for soap bubbles
p
H = − 2σ Laplace law for soap bubbles (4.54)
σ dA = −p dV
4π 3
σ d(4πR2 ) = −p d R
3
4π
σ 4π 2R dR = −p 3R2 dR
3
1 p
⇒ = − (4.55)
R 2σ
As the only compact surface with constant mean curvature (CMC-surface) is the
sphere, a soap bubble is spherical. CMC-surfaces are critical points of the area
functional under the constraint of a constant volume.
For p = 0 the shape equation is simply H = 0 , which describes a minimal
surface, i.e. a surface under tension with minimal energy given a particular
boundary curve. Those surfaces are always saddle-like, because H = 0 means
R1 = −R2 ⇒ K = − R12 , which is always negative. The implication is that those
1
surfaces cannot be compact, because a surface that is saddle-shaped cannot be
enclosed by a boundary. A well-known example for a minimal surface is the
catenoid connecting two wireframes, compare figure 4.12.
The solutions to the Euler-Lagrange equation for σ = 0 ´and finite κ are called
Wilmore surfaces. Because they are solutions to H = 2κ dA H 2 , minimal sur-
faces with H = 0 are included. But due to their saddle-like shape, minimal
4
ZC Ou-Yang and W Helfrich, Bending energy of vesicle membranes: General expressions
for the first, second, and third variation of the shape energy and applications to spheres and
cylinders, Phys. Rev. A 39: 5280-5288, 1989.
87
Figure 4.13: The shape diagram represents the minimal energy shapes for given points
in the phase space.
88
on the inner side of the vesicle, Nout is the number on the outside, and a, which
is the typical area per lipid, has the dimensions of nm2 . Since Ain and Aout do
not change, for energy reasons, ∆A0 stays constant for a vesicle. The bending
Hamiltonian for the ADE model is
α
ˆ
H = 2κ dA H 2 + (∆A − ∆A0 )2
2
The´ differential geometry result for the integrated mean curvature is ∆A =
2d dA H, with d = 4 nm. The resulting shape diagram is now two-dimensional
as shown in figure 4.13. It now contains the budded shape as well as non-
axisymmetric shapes like the starfish vesicle. In the literature, many similar
models have been discussed, including the spontaneous curvature and the bilayer-
couple models, to explain the zoo of vesicle shapes, but the ADE-model seems
to be the most appropriate one. Therefore it is also used as the standard start-
ing point to explain the shape of red blood cells, which are known to have very
asymmetric membrane leaflets.
89
Figure 4.14: Two dimensional shape diagram from ADE model. For each region, minimal
energy shapes are indicated. The horizontal axis is the reduced volume v; spheres again
correspond to v = 1. The vertical axis shows the effective differential area between inside
and outside monolayers.
90
Figure 4.15: Shape of red blood cells: comparison of experimentally observed and com-
putationally predicted shapes and their transitions. From the review by Gerald Lim and
colleagues.
91
therefore of oxygen in our body. During their 120 days lifetime, they travel 105
times through our circulation (each round trip takes 100 s) before they are sorted
out because they become stiffer. There are around 2.6 1013 RBCs in our body
(out of 3.1 1013 all together), making them the most abundant cell type8 . An
amazing number of 2 106 new ones are produced in every second in our bone
marrow. A RBC has a diameter of 8 µm, a thickness of 2 µm at the rim and of 1
µm at the middle of the biconcave disc. Its volume is 100 µm3 and its area 140
µm2 . This corresponds to a reduced volume of v = 0.642, in agreement with the
range from the vesicle theory in which we expect discocytes.
Under physiological conditions, area and volume do not change much and there-
fore can be taken as constant for our mathematical treatment. For area, this
results from the large area expansion modulus of KA = 0.5 J/m2 . The corre-
sponding energy is (KA /2)∆A2 /A0 and if we equate this with the bending energy
κ = 50 kB T of RBCs, we get ∆A/A0 = 10−5 , thus area does not change signifi-
cantly. In fact the membrane would rupture at one percent relative area dilation
and the large area expansion modulus protects it from this.
Volume control is more complicated. It mainly arises from osmotic pressure
arising from c0 = 290 mM of osmotically active molecules inside the cell. This
leads to an osmotic modulus KV = RT c0 = 7 105 J/m3 . Equating the energy
(KV /2)∆V 2 /V0 with the bending energy κ, we now get ∆V /V0 = 10−5 , thus
volume is also constant for practical purposes.
The standard model for RBC-shape was established in the beautiful paper by
Lim, Wortis and Mukhopadhyay in PNAS 2002. As shown in Fig. 4.16, the
plasma membrane of the RBC is reinforced by a polymer network (made mainly
from the intermediate filament spectrin) underlying it, thus forming a composite
or sandwich structure. The overall thickness however is so small that the system
can still be considered to be two-dimensional on the scale of the cell. Therefore the
authors expanded the ADE-model for the membrane by an elastic energy for the
polymer network. This elastic component is modeled as an isotropic hyperelastic
material. Isotropy is justified by the hexagonal network structure, but linearity is
not because the RBC is known to strain harden under the conditions in the blood
flow. Similar to the derivation of the Helfrich Hamiltonian, we write the elastic
Hamiltonian as a Taylor expansion, but this time not as a function of curvature,
but as a function of the two in-plane strain invariants α and β:
Kα
ˆ ˆ
H= dA α2 + α3 α3 + α4 α4 + µ dA β + b1 αβ + b2 β 2 (4.56)
2
where Kα is the stretch modulus and µ the shear modulus. The two strain
invariants follow from the principal extension ratios λ1 and λ2 of a deformed
ellipse as
1 λ 1 λ2
α = λ1 λ2 − 1, β = + −2 (4.57)
2 λ 2 λ1
In contrast to the Hamiltonian for the lipid bilayer, one now also needs a refer-
ence shape to calculate the elastic energy. A computational procedure has been
8
Compare the book by Ron Milo and Rob Phillips, Cell biology by the numbers, Garland
Science 2016
92
Figure 4.16: The shape of RBCs is determined by the nature of its composite membrane.
While the outside layer is a lipid membrane with differential lipid composition in the two
leaflets, the inside layer is a polymer network (made mainly from the polymer spectrin)
that is attached to the membrane at discrete points (band 3 tetramer, band 4.1 protein).
These anchor points form a hexagonal lattice and have a typical distance of 76 nm.
93
Figure 4.17:√Fluctuating membrane of lateral length L and typical deviation from a flat
membrane < h2 >
In this section we will investigate the mean square deviation < h2 > of a flat
lipid membrane fluctuating at temperature T, see figure 4.17. Its square root is a
measure for the size of typical excursions. In lowest order, the energy functional
for the almost flat membrane is (compare equations 4.26 and 4.27)
κ
ˆ ˆ
2
H [h(x, y)] = 2κ dA H = dxdy(∆h(x, y))2
2
1
ˆ
h(~x) = d~k h(~k) exp(i~k · ~x) (4.58)
(2π)d/2
1
ˆ
~
h(k) = d~x h(~x) exp(−i~k · ~x) (4.59)
(2π)d/2
ˆ∞
1
δ(k) = d~x exp(i~k · ~x) (4.60)
(2π)d
−∞
94
h(~x) has to be real
1
ˆ
h(~x) = d~k h(~k) exp(i~k · ~x)
(2π)d/2
1
ˆ
∗
= h(~x) = d~k h(~k)∗ exp(−i~k · ~x)
(2π)d/2
⇒ h(~k) = h(−~k)∗
The result is the same for k > 0 and for k < 0 and the case k = 0 is irrelevant,
because a(~k) = a(−~k) and b(~k) = b(−~k). Therefore we restrict the integration
to positive k, which gives a factor of 2. The bending energy is the sum of the
squares of the decoupled amplitudes. The k 4 -dependency is typical for bending.
The partition sum is a functional integral over all possible membrane conforma-
tions
ˆ
Z = Dh exp(−βH[h(x)])
Y
ˆ∞ ˆ∞
= da(~k) db(~k) exp(−βκk 4 [a(~k)2 + b(~k)2 ])
k>0−∞ −∞
Y kB T
=
k>0
κk 4
because this is simply a product of many Gauss integrals. For the free energy,
we therefore get
κk 4
ˆ
F = −kB T ln Z = kB T dk ln (4.62)
k>0 kB T
95
However, here we are interested in the correlation functions, not in Z or F directly.
For each ~k, there are two independent and harmonic degrees of freedom. We
therefore have
kB T
< a2 (~k) > =
2κk 4
kB T
< b2 (~k) > =
2κk 4
< a(~k)a(~k ′ ) > = < a(~k) >< a(~k ′ ) >= 0 etc.
~k 6= ~k′ , ~k 6= −~k′
Now the space-dependance drops out due to the delta function (the underlying
reason is translational invariance) and we are left with one integral only. If we
define a as the microscopic cutoff (molecular size) and L as macroscopic cutoff
(system size) and use d = 2, we get in polar coordinates:
2π
ˆa
2π kB T
< h2 (x, y) > = k dk
(2π)2 κk 4
2π
L
" 2 2 #
2π 1 kB T L a
= −
(2π)2 2 κ 2π 2π
= kB T
16π 3 κ
L2 =< h2 > (4.65)
96
Equation 4.65 shows that the mean square deviation is proportional to tempera-
ture T , inversely proportional to bending rigidity κ, and increases quadratically
with the system size L. Note that the limit a −→ 0 is unproblematic.
In order to better understand the fluctuations of membranes we can put in num-
bers:
κ = 20kB T
√
L = 10 nm ⇒ < h2 > = 1 Å
√
L = 1 cm ⇒ < h2 > = 100 µm
Thus the effect is relatively weak on the scale of vesicles, but relatively strong
on macroscopic scales. For a biomembrane fluctuations are relevant, but not on
small scales.
It is instructive to compare this result to the one for interfaces under tension (e.g.
oil droplets or soap bubbles). The we start from the Hamiltonian
σ
ˆ
H [h(x, y)] = dxdy(∇h(x, y))2
2
and therefore arrive at
kB T ~ ~ ′
< h(~k)h(~k ′ ) >= δ(k + k ) (4.66)
σk 2
The backtransform then gives
kB T L
< h2 >= ln( ) (4.67)
2πσ a
which is a much weaker dependance on L than for membranes. For σ = 100 erg/cm2
and a = 3 Å, a system size of L = 10 nm gives a mean deviation of 1.5 Å. For
L = 1 cm, this goes up only to 7.5 Å.
Another way to quantify membrane fluctuations is to investigate how much the
membrane looses its orientation due to fluctuations. A measure for this is the
persistence length Lp . Systems with characteristic length shorter than Lp can be
considered as elastic planes or rods (for polymers). The properties of systems with
characteristic length larger than Lp solutions can be described with statistical
methods for random walks. Formally, the persistence length is defined as the
length over which correlations in the direction of the normal are lost.
Let us again work in the Monge representation, see figure 4.18. The normal vector
is
1 hx 2
~n = √ hy , with detg = 1 + (∇h)
detg 1
1 hx
=q hy
1 + ∂x h2 + ∂y h2 1
Let us define the angle between normal vectors at different points on the mem-
brane as θ and
θ2 1 1
cos θ ⋍ 1 − = nz = q ≃ 1 + (∂x h2 + ∂y h2 )
2 1 + ∂x h2 + ∂y h2 2 | {z }
=(∇h)
2
97
Figure 4.18: Persistent length Lp for fluctuating membrane, Monge representation. ~n is
the normal to the membrane vector.
If we now set < θ2 >= π 2 for the extreme case that orientation has turned
around, we can define a length scale at which the membrane is not flat anymore:
Persistence length
Lp = a exp(2π 3 kBκT ) (4.68)
for membranes
The persistence length for membranes was calculated by de Gennes and Taupin
in 1982 10 .
For better illustration we look again at typical numbers. As already mentioned
in this section, for biomembranes κ ⋍ 20kB T , which makes Lp ≃ a · exp(400).
Although membranes are only 4 nm thick, this thickness is sufficient to conserve
their rigidity and flatness. Another example is the water-oil interface stabilized
by tensides (substances, that reduce surface tension and allow easier dispersion),
see figure 4.19. In this case κ ≃ 1 · kB T , Lp is small and the interface is thermally
roughened.
98
Figure 4.19: Surface interaction between oil and water mediated by tensides. The ten-
sides, represented by red circles with black tails, are responsible for the dispersion of oil
droplets into water. They reduce the surface tension, that is why the interface is rough.
(a) (b)
Figure 4.20: Steric interactions for a) stack of membranes and b) membrane trapped
between two walls. Characteristic dimension is the distance d for both cases. The
principal idea for the description of both cases is the same, but there are more degrees
of freedom for the stack of membranes.
99
density. We now will see that entropic effects lead to an effective repulsion be-
tween membranes. Consider a stack of membranes or a single membrane trapped
between two walls, see figure 4.20. We are interested in the free energy of the sys-
tem, which in this case is a function of the distance d between the membranes or
the membrane and the wall. We already can sense that this free energy increase
with increasing d because the membrane will gain entropy if the confinement
decreases.
Scaling argument
The membrane has “humps” of size h2 ∼ kBκT · L2p as calculated above. For each
hump, the membrane loses entropy per area kLB2 and the bending energy per area
p
2
h
for a hump is κ L2p
. From this argument we can conclude, that the free energy
per area is
!2 !
∆F h −kB (kB T )2 1
∼κ −T ∼ · 2
A L2p L2p κ h
For a membrane in a stack or between two walls, h scales like d and therefore we
get the fluctuation or Helfrich interaction11 :
∆F (kB T )2 1
A ∼ κ · h2
Helfrich 1978
Although this argument involves bending energy, this too arises from thermal
fluctuations. Therefore the whole effect is a thermal one and vanishes with T → 0.
An exact solution is not known, but a reasonable calculation starts with the
confinement effect modeled by a harmonic potential12 . Thus we consider the
mean squared deviation in the Monge representation for an almost flat membrane
that fluctuates under a harmonic potential:
κ 1 2 1
ˆ ˆ n o
2
H= dx dy (△h) + 4 h = dx dy κ(△h)2 + γh2
2 ξ 2
1 kB T
ˆ
2
<h > = 2
2π dk k 4
(2π) κ(k + ξ −4 )
kB T 2 kB T
= ·ξ = √
8κ 8 κγ
11
Helfrich, W. Steric interaction of fluid membranes in multilayer systems. Z. Naturforsch 33
(1978): 305-315.
12
Janke, W., and H. Kleinert. Fluctuation pressure of membrane between walls. Physics
Letters A 117.7 (1986): 353-357.
100
We now assume a simple geometrical scaling of the excursion with the confinment,
< h2 >= µd2 . Here µ is a constant prefactor, that has been found in Monte-Carlo
computer simulations to be µ = 16 . Combining the two expressions for < h2 >,
we can solve for ξ as a function of d. Because we have a harmonic (Gaussian)
system, for the free energy difference between the confined and the free membrane
we get (compare the introduction, free energy of a harmonic system)
∆F
= −kB T · ln Z
A !
kB T k 4 + ξ −4
ˆ
= 2π k dk ln
(2π)2 k4
kB T −2
= ·ξ
8
From this follows:
∆F (kB T )2
A = 64κµd2
Steric interaction between membranes
This is the same result as from the scaling analysis, but now with exact pref-
actors. This result has been confirmed both with computer simulations and in
experiments.
If we repeat the same analysis for the case of surface tension, we have for the
squared mean displacement
1 kB T
ˆ
2
< h >= 2
2π dk k 2 (4.69)
(2π) σk + γ
Now the integral is not (1/2) arctan(k 2 ), but (1/2) ln(1 + k 2 ), thus it diverges
for large k and we have to use a microscopic cutoff. If we combine both surface
tension and bending rigidity, however, we get a well-defined result again.
kB T
< h(~k)h(~k ′ ) >= δ(~k + ~k ′ ) (4.70)
σk 2 + κk 4 + γ
13
Nir S Gov and Sam A Safran, Red blood cell membrane fluctuations and shape controlled
by ATP-induced cytoskeletal defects, Biophysical journal 88.3 (2005): 1859-1874.
101
This procedure has been applied successfully to RBCs under various conditions
and is has been found that shape is the main determinant of the fluctuations14 .
We note that assuming an almost flat membrane is a strong assumption and that
a more rigorous analysis had to consider also the role of curvature.
Alternatively, one can assume that the whole shell is one composite and fluctuates
as such, as we have assumed above to derive the minimal energy shape. Then one
has to work with thin shell elasticity and the results are much more complicated.
In this way, it has been shown that at low and high frequencies, the fluctuations
are dominated by active and passive fluctuations15 . Active fluctuations depend on
ATP and arise e.g. from the actin-spectrin network or ion pumps and channels.
14
Yoon, Young-Zoon, et al., Flickering analysis of erythrocyte mechanical properties: depen-
dence on oxygenation level, cell shape, and hydration level, Biophysical journal 97.6 (2009):
1606-1615.
15
Turlier, Herve, et al., Equilibrium physics breakdown reveals the active nature of red blood
cell membrane fluctuations, Nat. Phys. (2016). There are two excellent reviews on this sub-
ject: Turlier, Herve, and Timo Betz, Fluctuations in active membranes, Physics of Biological
Membranes, Springer 2018. 581-619, and Turlier, Herve, and Timo Betz, Unveiling the Active
Nature of Living-Membrane Fluctuations and Mechanics, Annual Review of Condensed Matter
Physics 10 (2019): 213-232.
102
Chapter 5
Physics of polymers
Polymers are chain molecules that can be described as space curves in three di-
mensions ~r(s) using the language and tools of differential geometry as introduced
in the membrane chapter. Motivated by the phenomenological approach to mem-
branes, we could start in a continuum framework with a bending Hamiltonian:
ˆL !2
κp d2~r(s)
H[~r(s)] = ds (5.1)
2 ds2
0
103
H H polymerization H H
C C C C
H H H H N
ethylene polyethylene
Figure 5.1: Polymerization of ethylene. Polyethylene (PE) is made by opening the double
bond between the carbon atoms in ethylene, flipping it over and thus connecting to the
next ethylene monomer. The subscript N is the degree of polymerization.
polyvenylchloride PVC
H H Cl
C C C C C C versus C C C C
H Cl Cl Cl Cl
venylchloride
(a) (b)
104
Type of polymer Sketch Example
Homopolymers are mostly
homopolymer -A-A-...-A- synthetic polymers, e.g. PE
DNA, which has 4 different monomers
heteropolymer -A-B-A-C-...- Proteins, which have 20 different monomers
Those are heteropolymers with two blocks,
diblock-copolymer -A-...-A-B-...-B- each build up of a different monomer.
This structure is similar to lipids.
H-polymer
ladder polymer
105
Figure 5.3: The elastic properties of polymer gels can be studied by putting the them
between two walls and then rotating or shearing those walls against each other. The
typical mechanical behaviour of the polymer network is depicted on the right side. The
logarithm of the elastic moduli is shown as a function of the logarithm of the frequency
w. At low frequencies the material is viscous, but it becomes elastic at high frequencies.
The R~ i are the position vectors for the nodes of the chain. R
~ is the end-to-end
vector, giving a characteristic dimension of the polymer, see figure 5.4. As each
106
Figure 5.4: Freely jointed chain FJC model for polymers. A short polymer is schemati-
cally depicted, as a chain consisting of segments r~i , represented as vectors. All segments
have the same length a. R ~ is the end-to-end vector.
N
X
< R >= < r~i >= 0 .
i=1
In analogy to the mean squared deviation < h2 > for membranes, we look at the
mean squared end-to-end distance
!
X X
< R2 > = < r~i r~j >
i j
N
X X
= < r~i 2 > + < r~i · r~j > = N a2 (5.4)
| {z } | {z }
i=1 i6=j
=a2 =0
√ typical extension of
R= Na polymer chain (5.5)
The square root relation is typical for a random walk. We introduce time t = N τ
(with stepping time τ ) and get
with D = a2 /2τ the diffusion constant and d spatial dimension. In fact our
polymer model is exactly the prescription of how to implement a random walk.
In a real polymer there are correlations between the different bond vectors, <
r~i · r~j >6= 0 even for i 6= j. However, in most cases, the polymer becomes “ideal”
in the sense that there are no correlations between monomers at large distance
along the chain, < r~i · r~j >= 0 for |i − j| −→ ∞. Therefore the sum over these
107
Ideal polymer C∞ b[Å]
polyethylen −CH2 CH2 − 7.4 14
polybutadiene −CH2 CH = CH CH2 − 5.3 9.6
polyisoprene (rubber) −CH2 CH = CH CH CH3 − 4.6 8.2
polydimethylsiloxane (elastomere) −O Si (CH3 )2 − 6.8 13
Table 5.3: Flory’s characteristic ratio and Kuhn lengths for different polymers
= CN N a2 N−→
−→∞
C∞ N a2 (5.6)
with C∞ = Ci ∀i with 1 ≤ C∞ < ∞. C∞ is called “ Flory’s characteristic ratio”,
see table 5.3.
Ideal polymers correspond to a FJC with redefined monomer length b and degree
of polymerization N :
L = N · b, < R2 >= N · b2 = b · L
<R2 >
b=
<R2 >
L
L2 Kuhn length (5.7)
N= b2
= <R 2>
The Kuhn length is a measure for the statistical segment length and tabulated
in table 5.3.
108
Figure 5.5: Schema of the freely rotating chain model (FRC). Here the length a and the
bond angle θ, between the segments, are kept constant. The torsion angle φ is still free
and makes the FRC flexible.
The persistence length has the same meaning as in membrane physics, it denotes
the length scale over which the correlations decay.
We now can use this exponential decay to calculate the mean squared end-to-end
distance:
N X
X N
< R2 > = < r~i · r~j > (5.10)
i=1 j=1
N
X i−1
X N
X
= < r~i · r~j > + < r~i >2 + < r~i · r~j > (5.11)
i=1 j=1 j=i+1
N
X i−1
X N
X
= a2 N + a2 (cos θ)i−j + (cos θ)j−i (5.12)
i=1 j=1 j=i+1
N i−1 N −i
!
X X X
2 2 k k
= a N +a (cos θ) + (cos θ) (5.13)
i=1 k=1 k=1
The two sums can be extended to infinity because at large distances, the corre-
lation has decayed. We then simply have a geometrical series:
∞
X cos θ
(cos θ)k =
k=1
1 − cos θ
Therefore
cos θ
< R2 >= a2 N + 2a2 N (5.14)
1 − cos θ
109
such that the final result reads
mean square
1+cos θ
< R2 >= N a2 1−cos θ
end-to-end distance (5.15)
for freely rotating chain
1+cos θ
If we compare this result with equation 5.6 , we see that C∞ = 1−cos θ . The
values for C∞ are typically between 5 and 8, see table 5.3.
110
(a) (b)
Figure 5.6: a.) The worm-like chain model describes an elastic rod. This polymer model
is similar to the Helfrich Hamiltonian for membranes. b.) The dependency of the mean
squared end-to-end distance on the ratio of the contour and persistence lengths of a
polymer in the WLC model.
The general expression is a smooth crossover between the two, see figure5.6b:
< R2 > 2 2
2
= − 2 (1 − e−x )
L x x
with x = L/lp . This defines three different regimes: the flexible chain at x ≫ 1,
the semiflexible chain with x ≈ 1, and the rigid polymer with x ≪ 1. Biological
examples are DNA, actin and microtubules.
The mean squared end-to-end distance < R2 > gives a measure for the extension
of the polymer, but it is very hard to measure it directly and what one usually
measures in experiments is the radius of gyration (e.g. with light scattering or
size-exclusion chromatography). We now clarify the relation between the two.
Assumed the monomer mass is constant, for the center of mass we have:
N
~ cm = 1
X
R ~i
R
N i=1
111
Now we will calculate the mean squared radius of gyration < Rg2 >:
N
1 X
Rg2 = ~ cm )2
~i − R
(R
N i=1
N
1 X ~ i 2 − 2R
~iR ~2 )
~ cm + R
= (R cm
N i=1
N N N N
1 X
~i2 1
X 1 X ~i 1
X
~j
= R − 2R R
N i=1
N j=1
N i=1
N j=1
| {z }
=1
N
! N
1 X
~i 1 X ~j
+ R R
N i=1
N j=1
1 XX ~ 2 ~iR
~j + R
~iR
~ j)
= (Ri − 2R (5.22)
N2 i j | {z }
~i R
=−R ~j
This expression does not depend on the choice of summation indices and we
rewrite it in a symmetric form:
1 1 X X ~ 2 ~ ~ XX
~j 2 − R
Rg2 = 2
(Ri − Ri Rj ) + (R ~ i )
~j R
N 2 i j j i
1 XX ~ 2 ~ jR ~j 2 )
~i + R
= (Ri − 2R
2N 2 i j
1 XX ~ ~j )2
= (Ri − R
2N 2 i j
N X N
1 X ~j )2
~i − R
= (R
N 2 i=1 j=i
1
< Rg2 >= ~j )2 >
PN PN ~i − R
N2 i=1 j=i < (R (5.23)
The radius of gyration can be expressed in terms of the average square distance
between all pairs of monomers.
For an ideal polymer chain, the sums can be changed into contour integrals:
ˆN ˆN
1
< Rg2 >= 2 du ~
dv < (R(u) ~
− R(v))2
>
N
0 u
We now use the fact that the contour between u and v should behave also like a
(shorter) ideal chain, see figure 5.7:
2
~
< (R(u) ~
− R(v)) >= (v − u)b2
112
Figure 5.7: Integration along the ideal polymer. Assumption that the contour between
u and v also behave like an ideal chain.
<R2 >
< Rg2 >= 6 Debye result (5.25)
This is a very important result. It means that for an ideal chain we can work
both with < R2 > or < Rg2 >, they are essentially the same, except for a constant
factor of 6.
113
with
N
1
X
ˆ P
~ ~ r~i )
~−
δ(R r~i ) = 3
d~k eik(R−
i=1
(2π)
We therefore obtain
ˆ N
1 1
ˆ
~~ −i~k~
~ =
p(R) d~k eikR d~r e r
δ(|~r| − a)
(2π)3 4πa2
ˆ∞ ˆ2π ˆ1
1 1
ˆ
−i~k~ 2
d~r e r
δ(|~r| − a) = r dr dφ d(cos θ) e−ikr cos θ δ(r − a)
4πa2 4πa2
ˆ0 0 0
2π sin kr
= 2
dr r2 δ(r − a)2
4πa kr
sin(ka)
=
ka
because
1
1 sin(kr)
ˆ
e−ikrx dx = (e−ikr − eikr ) =
−1 −ikr kr
For N ≫ 1 we can assume N −→ ∞, so only ka ≪ 1 is relevant. We now make
a Taylor expansion around 0:
!N
k 2 a2
N
sin(ka) N k2 a2
≈ 1− ≈ e− 6
ka 6
So we get
1
ˆ 2 2
~~ Nk a
~ =
p(R) d~k eikR e− 6
(2π)3
| {z }
Gauss integral
1
ˆ
Y 2 a2
= dkα e−ikα Rα −N kα 6
(2π)3 α=x,y,z
3
1 6π 2 ~2
3 R
= e− 2 a2 N (5.28)
(2π)3 N a2
The distribution function of the end-to-end vector is Gaussian. The feature that
R > N a is an artifact of our expansion, but that does not matter, because the
respective weights are negligible.
114
(a) (b)
Figure 5.8: (a) Gaussian distribution of one component Rα of the end-to-end distance
~ (b) Probability distribution of the radial component of R
vector R. ~ (in spherical coordi-
nates). Note the similarity to the Maxwell-Boltzmann distribution.
Figure 5.9: The Gaussian Chain Model. The Gaussian length distribution of the bond
lengths are depicted as springs with the entropic spring constant k.
In Cartesian coordinates, equation 5.29 reads for the single components (compare
also figure 5.8a):
!1 !
~ = Q 3 2
3Rα2
p(R) α=x,y,z exp − (5.30)
2πN a2 2N a2
dRα p(Rα ) = 1
´
⇒
N a2
< Rα2 >= dRα p(Rα )Rα2 =
´
⇒
3
In spherical coordinates, one finds for the modulus of the radius:
!3 !
3 2
3R2
p(R) = exp − · 4πR2 (5.31)
2πna2 2N a2
This result for p(R) (figure 5.8b) is equivalent to the Maxwell-Boltzmann distri-
bution for the distribution of the modulus of the velocity for an ideal gas. The
same Gaussian distribution is obtained by starting from a symmetric random
walk on a lattice, i.e. from the binomial distribution.
These results from the FJC motivate us to define a new polymer model that as-
sumes the Gaussian property to be valid for every segment. In the Gaussian chain
115
model (GCM) one assumes that every bond has a Gaussian length distribution
(instead of a fixed length a as in the FJC):
!
3~r2
3
3 2
p(~r) = exp − 2 (5.32)
2πa2 2a
with
where we have used the substitution ds = adn. Note that this Hamiltonian is
fundamentally different from the WLC Hamiltonian from equation 5.1 because it
describes stretching and not bending.
We now consider the free energy of the Gaussian chain:
F = |{z}
U − TS
=0
~
= −T · kB ln Ω(R)
ˆ
~ ·
= −kB T · ln p(R) ~ Ω(R)
dR ~
Since ~ Ω(R)
~ is independent of R,
~ the free energy F can be written as:
´
dR
3 R~2
F = kB T + F0 (5.36)
2 N a2
where F0 does not depend on R. ~ The free energy of a Gaussian chain increases
quadratically with R, because the number of possible configurations and hence
the entropy decreases. This leads to Hooke’s law:
3kB T ~
F~ = − ·R (5.37)
N a2
116
Figure 5.10: A force applied to two beads attached to the polymer, e.g. by optical
tweezers. The polymer is stretched and the force needed to stretch it to a certain length
is measured.
where 3kB T /(N a2 ) is the entropic spring constant of the whole chain.
Note that for higher temperatures, the entropic spring constant increases or, in
other words, the bonds become harder to stretch. In the limit T → ∞, the chain
contracts into a single point. The reason for this suprising behaviour is that the
effective energy that is needed to stretch the polymer is entirely related to the
loss of entropy. It is therefore easier to stretch polymers with a larger number of
monomers N , larger monomer size a and lower temperature T . This is different
for energy-dominated materials such as metals, which become softer at higher
temperature.
The entropic spring constant energy of the Gaussian Chain suggests to study
the behavior of a FJC under stretch. Imagine placing beads at the ends and
pulling them apart along the z-axis with optical tweezers (figure 5.10). Today,
the pulling of biopolymers with AFM, optical or magnetic tweezers, electric fields
or hydrodynamic flow, to name but a few, is a standard experiment in biophysics.
Obviously, the Gaussian result (equation 5.29) cannot be true for large extensions,
i.e. close to the contour length. In the following we will approach the problem
of calculating the force-extension curve for finite contour length first by a scaling
argument and then by an analytical calculation.
Scaling analysis
ζ 2 = gb2 (5.38)
where ζ denotes the blob size and g denotes the number of monomers per blob.
Hence, the total number of blobs is simply N/g. From here on we use the symbol
b for the Kuhn length as the segment length.
117
Figure 5.11: Polymer depicted as a chain of "blobs" which on the blob scale ζ behaves as
an unperturbed random walk. On the large scale Rz , the blobs are oriented in z-direction.
118
Figure 5.12: A force Fz applied to a freely jointed chain. The chain consists of N bond
vectors ~ri with a fixed length b.
N ˆ 2π ˆ 1
Y − k HT
⇒Z = dφi d(cos Θi ) e B
| {z }
i=1 0 −1
:=xi
N ˆ
Y 2π ˆ 1 F b X N
z
= d(xi ) exp
dφi · xi
kB T
i=1 0 −1
| {z } i=1
:=f
ˆ 1
= (2π dx ef x )N
−1
2π f 4π
= [ (e − e−f )]N = ( sinh(f ))N (5.45)
f f
With the free energy, the expectation value of the spatial extension in z-direction
can be computed:
∂F ∂f
< Rz >= − ∂F z
= − ∂F
∂f · ∂Fz
where we introduced the Langevin function L(f ). Equation 5.47 has two inter-
esting limits (compare also figure 5.13):
119
In this regime L(f ) can be approximated by a linear function:
ef + e−f
coth f =
ef − e−f
1 f
(1 + f + 21 f 2 ) + (1 − f + 12 f 2 ) f + 2
≈ =
(1 + f + 21 f 2 + 16 f 3 ) − (1 − f + 21 f 2 − 61 f 3 ) 1 + 16 f 2
1 f 1
≈ ( + )(1 − f 2 )
f 2 6
1 1
≈ + f
f 3
f
⇒ L(f ) ≈
3
Fz b
⇒ < Rz >= bN (5.48)
3kB T
ef
coth f =1
≈
ef
1
⇒ L(f ) = 1 −
f
!
kB T
⇒ < Rz >= bN 1 − (5.49)
Fz b
where ~t is the tangential vector (|~t| = 1). There are three ways to derive this
Hamiltonian:
120
Figure 5.13: The expectation value of the extension in z-direction as a function of the
(dimensionless) force parameter f = Fz b/(kB T ) (left) and vice versa. For small force
and small extension, respectively, the relation is approximately linear. For large forces
the extension approaches the contour length L and the force diverges.
2. From beam theory: With the Young’s modulus E of the elastic rod one
finds
κp = E · I (5.52)
where
r2 R
r2 πR4
ˆ ˆ
I = dA = 2π rdr = (5.53)
2 0 2 4
The bending energy can be calculated by discretizing the WLC into small seg-
ments as shown in figure 5.14. This model is similar to the FRC. However, the
main difference is that for the FRC the angle Θ was held constant whereas now
the second moment < Θ2 >∼ kB T b/κp is a constant. The bending Hamiltonian
121
in equation 5.51 then gives
N
X κp
Eb = · (~ti − ~ti−1 )2
i=1
2b
N
X κp
= · (1 − cos Θi )
i=1
b
N
small curvature X κp 2
⇒ Eb ≈ Θ (5.54)
i=1
2b i
Persistence length
where
κp
lp = persistence length (5.56)
kB T
with f (x) = 2x − 2x2 (1 − e−1/x ) and x = lp /L. Below we will make us of the
scaling function f (x).
122
Extension in z-direction
1 d ln Z
ˆ H
−
⇒ extension < Rz > = D~t Rz e kB T = (5.59)
Z df
In contrast to the FJC an exact solution to equation 5.59 is not known. However,
the two asymptotic limits can be treated analytically:
Since
L L
1 ~ 2 > ≈ p 1 · 2Llp
ˆ ˆ
L≫l
du dv < tz (u)tz (v) >0 = <R
0 0 3 3
Therefore the extension of the WLC in response to a small stretch exhibits, similar
to the FJC, a linear force-extension dependency with an entropic spring constant
kW LC = 3kB T /(2lp L). Recall, that for the entropic spring constant for the FJC
we found kF JC = 3kB T /(bL) (compare equation 5.48).
123
In this regime we are dealing with an almost straight chain and can therefore use
a Monge parametrization for ~t (compare figure 5.15a):
tx
~t =
ty
(5.62)
1 2 2
1 − 2 (tx + ty )
where for the z-component we have used the fact that ~t is normalized and a
Taylor expansion. Our Hamiltonian now reads
!2 !2 !2
H lp dtx dty dtz
ˆ ˆ
= ds [ + + ] − f ds tz
kB T 2 ds ds ds
| {z }
≈0
" 2 2 #
lp dtx dty f
ˆ ˆ h i
= ds + + ds t2x + t2y − f · L (5.63)
2 ts ds 2
This Hamiltonian is quadratic and therefore the partition sum is a Gaussian path
integral (the constant term does not matter). In Fourier space we thus have
Fourier kB T
⇒ |tα (k)|2 = (5.64)
kB T (lP k 2 + f )
L ˆ L
1
ˆ
⇒< Rz > = ds < tz >= ds < (1 − (t2x + t2y ) >
0 0 2
ˆ L
1 1
= L− ds (< t2x > + < t2y >) = L − L · 2 < t2x >
2 0 2
ˆ ∞ ˆ ∞
1 1 1 1
= L · (1 − dk 2
) = L · (1 − dk q )
2π −∞ lp k + f 2πf −∞ l
( p k)2 + 1 f
s
∞
1 f 1
ˆ
= L · (1 − dk′ )
2πf lp −∞ k′2 +1
| {z }
π
1
= L(1 − p ) (5.65)
2 f lp
⇒ L−<Rz >
L = √1 (5.66)
2 f lp
√
This is a square-root divergence ∼ 1/ Fz (figure 5.15b). Recall that the FJC
in the large extension regime diverges with 1/Fz (compare equation 5.50) which
is crucially different. In experiments, stretching semiflexible biopolymers like
dsDNA has shown that they can not be described by the FJC2 .
2
Smith, Steven B., Laura Finzi, and Carlos Bustamante. Direct mechanical measurements
of the elasticity of single DNA molecules by using magnetic beads. Science 258.5085 (1992):
1122-1126; Bustamante, C., Marko, J. F., Siggia, E. D., and Smith, S. Entropic elasticity of
lambda-phage DNA. Science (1994): 1599-1599; Bustamante, Carlos, Zev Bryant, and Steven
B. Smith. Ten years of tension: single-molecule DNA mechanics. Nature 421.6921 (2003):
423-427.
124
(a) (b)
Figure 5.15: a.) A WLC under large stretch. The polymer is now almost straight and can
be assumed to have no overhangs. Hence, in analogy to membrane physics, we can chose
a parametrization that is similar to the Monge parametrization. b.) Force-extension
dependence for the FJC and the WLC. For the FJC the scaling is (L− < Rz >)/L ∼ Fz−1 ,
whereas for the WLC we find (L− < Rz >)/L ∼ Fz .
−1/2
Although an exact formula for the WLC is still lacking, the two limits shown here
can be combined in an interpolation formula with an error smaller than 10%3 :
3 < Rz >
1. small stretch f · lp =
2L
1
2. large stretch f · lp =
4 · (1 − <RLz > )2
The two limiting cases of the stretched WLC can be obtained also from a blob
scaling analysis, which helps to better understand the underlying physics. We
consider a chain segment of length l which is bent to an angle θ. As discussed
above, this costs the bending energy
κp (1 − cos θ) κp θ 2
Eb ∼ ≈ (5.68)
l l
which diverges with l → 0, because we would get infinite curvature. In order to
work against the external force Fz , we need the stretching energy
that increases with l. Therefore a crossover length ξ exists at which the two
energies balance: s
r
κp kB T lp
ξ := = . (5.70)
Fz Fz
3
JF Marko and ED Siggia: "Stretching DNA", Macromolecules 1995, 28:8759–8770
125
We interpret ξ as the contour length per blob. Below it, the chain does not feel
the effect of force and is dominated by bending. Above it, the chain becomes
elongated in z-direction and is dominated by stretching. In the blob picture, we
assume that we have an unperturbed WLC below ξ and a stretched FJC of blobs
above ξ.
We next recall the two scaling functions that we have calculated above. For the
unperturbed WLC, we have defined a scaling function f (x) for the mean squared
end-to-end distance in eq. 5.57. We now use it to define the size of a blob as
lp
b2b := ξ 2 f ( ) . (5.71)
ξ
For the FJC of blobs, we can use the Langevin function L(x) defined in eq 5.47.
We note that we have L/ξ blobs, each of size bb , and therefore the overall relative
extension will be
< Rz > 1L Fz b b
∼ bb L( ) (5.72)
L Lξ kB T
where bb depends on the regime in which the scaling function f (x) is evaluated. A
closer look shows that the overall result is controled only by one scaling parameter,
namely f := Fz lb /(kB T ).
We now can look at the two limiting cases. For strong stretching, f ≫ 1, we have
L(x) = 1 − 1/x and f (x) = 1. Thus bb = ξ (the blob is rigid with linear scaling)
and therefore
< Rz > kB T
∼1− . (5.73)
L Fz ξ
We rearrange to find s
L− < Rz > kB T kB T
∼ (5.74)
L Fz ξ F z l p
exactly as found above, except that the scaling analysis misses a factor of 2.
For weak stretching, f ≪ 1, we have L(x) = x/3 and f (x) = 2x. Thus b2b = 2ξlp
(the blob is flexibel with square root scaling) and we get
p p
< Rz > ξlp Fz ξlp Fz l p
∼ = (5.75)
L ξ kB T kB T
because ξ cancels out. This is the linear response regime that we also found
above. Here we miss a numerical factor of 2/3 compared with the exact result.
Overall we conclude that the blob analysis gives the right scaling results in both
limits, in particular the inverse square root for the divergence at strong stretching,
which sets the WLC apart from the FJC, and the linear response regime at weak
stretching.
126
but also for other biopolymers like actin, it is known that after the thermal
fluctuations in the contour length have been pulled out, the backbone can give out
additional length due to internal changes (overstretching in the case of dsDNA,
twist in the case of actin). This situation is described by the stretchable WLC-
model, which can be solved with the same methods as described above for the
WLC-model, and which is a combination of the GC and the WLC.
The following references are recommended for further reading:
• R Phillips et al., Physical biology of the cell, chapter 10; especially appendix
10.8 on the math of the WLC on page 401
• P Nelson, Biological Physics, very detailed discussion of different models
• Kroy, Klaus, and Erwin Frey. Force-extension relation and plateau modulus
for wormlike chains. Physical Review Letters 77.2 (1996): 306.
• J Kierfeld et al. Stretching of semiflexible polymers with elastic bonds, Eur.
Phys. J. E 2004, 14:17-34
• Koester, S., J. Kierfeld, and T. Pfohl. Characterization of single semiflexible
filaments under geometric constraints. The European Physical Journal E
25.4 (2008): 439-449.
127
Figure 5.16: The extension of DNA has been measured on supported bilayers and resulted
in an exponent 0.79 very close to 2D Flory theory. From B. Maier and J.O. Rädler, Phys.
Rev. Lett. 82, 1911, 1999.
128
Figure 5.17: Experimental results and scaling laws for the modulus of different polymer
networks. (a) Note that the synthetic polyacrylamide gel is the only one that does not
strain-stiffen. (b) Actin network crosslinked by scruin. (c) Neurofilament network. (d)
Polyisocyanopeptide hydrogel. The power 3/2 is the prediction of the affine thermal
model. Taken from Broedersz and MacKintosh review, figure 15.
The mechanical stability of cells and tissues results mainly from networks of
semiflexible polymers (e.g. actin inside the cells and collagen between the cells).
These kinds of networks are stabilized both by topological entanglement and by
crosslinkers (e.g. alpha-actinin, fascin, filamin, fimbrin, scruin etc for actin).
Despite the fact that molecular details and network architecture can vary widely
in these systems, they all share one outstanding property, namely that they stiffen
under strain, as shown in the experimental plots shown in figure 5.17. We have
seen this already for the single WLC, but it is non-trivial to find this result also
for the network. A nice review on this subject is by Chase P. Broedersz and Fred
C. MacKintosh, Reviews of Modern Physics 2014, volume 86, pages 995-1036.
We cannot go into the details here, but would like to mention one difficulty: if
one couples different polymers into a bulk material, most deformation modes will
include both stretch and compression of polymers. However, these two modes are
very asymmetric on the level of the molecules. Under compression, the polymer
does buckle at the Euler threshold. This can be seen easily by noting that the
thermal fluctuation of a beam is
H ∼ (σ + κk 2 )k 2 (5.82)
which can become negative for σ < −κk 2 . Because the critical wavelength is
129
related to beam length L by k ∗ = π/L, we get for the critical tension at buckling
π R
σc = κ( )2 = π 2 E( )2 (5.83)
L L
with bending stiffness κ = ER2 , Young’s modulus E and radius R. Thus the
longer a polymer, the more easily it buckles. A complete theory of a poly-
mer gel has to incorporate this asymmetry, the scale on which the polymers
are crosslinked, and the nature of the crosslinks.
130
Chapter 6
Molecular motors
Molecular motors are molecules that generate motion and force. They do this
by converting electrochemical energy into mechanical work, for example by hy-
drolysing ATP or by letting ions flow down a gradient. Thus they work like heat
engines, but they cannot be Carnot engines, because molecular diffusion is too
fast as to allow for any temperature gradients. Thus they have to achieve the
conversion without the intermediate form of heat and to operate at constant tem-
perature (isothermally). Molecular motors are extremely fascinating molecular
machines and it is still not completely clear how they have been optimized by
evolution to perform their tasks. An important aspect of understanding them
is to build new ones, for example by reengineering their different parts or by
using different material (e.g. small molecules or DNA rather than proteins). In
2016, the Nobel prize for chemistry has been awarded for the design and synthe-
sis of molecular machines and a Nobel prize for medicine for the investigation of
biological molecular motors might come in the future.
Why did nature evolve motors ? Obviously this is a very direct way to generate
force, e.g. in the muscle for moving body parts or in the beating flagella of sperm
cells. In regard to transport, for examples of vesicles and organelles, but also
of viruses, motors are needed not only to provide specificity and direction, but
also to beat the physical limits of diffusion. With a typical molecular diffusion
constant of (10 µm)2 /s, diffusion starts to become slow in regard to the required
response times of s on the length scale of cells (10 µm). With a typical velocity of
µm/s, molecular motors outcompete diffusion on cellular and tissue length scales.
However, we also note that for body length scales, we need other transport modes.
For example, small molecules such as hormones and many cell types (red blood
cells, platelets, white blood cells, stem cells and cancer cells) are transported with
the blood (average velocity 0.4 m/s in the aorta and 0.3 mm/s in the capillaries)
and nerve impulses are transmitted as action potentials (velocities 10-100 m/s).
In this chapter, we will discuss the theoretical basis of understanding molecu-
lar motors. As we will see, the appropriate framework is the one of stochastic
equations (master equation, Fokker-Planck equation) and the theory of molecular
motors has advanced considerably over the last two decades and still is a very
131
active research area 1 .
6.1 Classification
Molecular motors can be classified in the following way:
Translational motors These motors move along tracks, e.g. myosin motors
along actin filaments (e.g. in the muscle), kinesin and dynein along micro-
tubules (e.g. kinesin in axons for transport towards the synaptic cleft, and
dynein in cilia and flagella to bend them), and polymerases and helicases
along DNA.
Rotary motors These motors typically have a stator embedded into the mem-
brane and containing a rotor. The most important example is the F0 F1
ATP Synthase, which in cells of all species generates ATP from ADP and
Pi (1 ATP per 120 degree rotation, at 100 Hz this gives 300 ATP per sec-
ond). It is driven by a proton gradient and needs six protons for each turn.
An adult burns 120 W and needs 2.400 kcal / day and thus 1.7 × 1026
ATP molecules, amouting to 140 kg that are essentially produced in our
mitochondria. The required energy comes from our metabolism (aerobic
respiration of glucose, which essentially was produced before by plants us-
ing photosynthesis). If there is plenty of ATP, the motor reverses and builds
up the proton gradient. Another famous example is the bacterial flagellar
motor, which is basically constructed like a ion turbine using 1.000 protons
to drive one turn. This motor has to create more torque than the ATP
Synthase because it has to turn the bulky flagellum.
Although these motors are very different on the molecular level, they share the
basic principle, namely stochastic operation at constant temperature to create
biased movement along cyclic phase space trajectories.
1
In the introduction and the discussion of the force-velocity relation, we follow the book
Physical biology of the cell by Rob Phillips and coworkers. For the more mathematical discus-
sion, we follow two excellent review papers on this subject: Frank Jülicher, Armand Ajdari and
Jacques Prost, Modeling molecular motors, Reviews of Modern Physics 69, 1269-1281, 1997;
Tom Duke, Modelling motor protein systems, course 3 in Physics of bio-molecules and cells,
volume 75 of the Les Houches series, pages 95-143, 2002.
132
In order to foster model building, we start with the simplest example, namely a
translational motor walking along a track. We consider a processive motor like
kinesin or myosin V, that can make many steps without falling off the track (this
is different for non-processive motors like myosin II, that stay on track only for
a short time and thus can work productively only in groups). Such motors are
typically two-headed and move in a hand-over-hand fashion. Moreover each step
is related to exactly one ATP being consumed. We label the track position by
the spatial coordinate x and assume that each motor has only a finite number of
discrete states, which we label with the index m. Thus our central quantity is
the probability pm (x, t) to be in state m and at position x at time t.
The two gain terms come from motors hopping in from left and right, respectively,
and the two loss terms come from motors hopping away to the left and right,
respectively. We rearrange and take the continuum limit ∆t → 0 to get
We next take the continuum limit in regard to space and use the Taylor expansion
1
p(x ± a, t) ≈ p(x, t) ± p′ (x, t)a + p′′ (x, t)a2 (6.3)
2
We then end up with the famous Fokker-Planck or Smoluchowski equation
a2
v = (k+ − k− )a, D = (k+ + k− ) (6.5)
2
For a Delta function as initial condition, this equation is solved by
1 2
p(x, t) = √ e−(x−vt) /4Dt (6.6)
4πDt
Thus the motor moves with a drift velocity v to the right, but it also disperses
with a diffusion constant D.
133
A B stall force
force
C D
Figure 6.1: Force dependence. (A) Scheme how force will change the free energy land-
scape of a motor hopping to the right. (B) Force-velocity relation when force dependence
is in k+ . (C) Force-velocity relation when force dependence is in k− . (D) Some exper-
imentally measured force-velocity relations: kinesin (green), RNA polymerase (blue),
phage packaging motor (red). All four graphs taken from the book Physical Biology of
the Cell, chapter 16 on molecular motors.
We also note that one can derive a dispersion relation from here. We use the
Fourier ansatz p(x, t) = C exp(i(kx − ωt)) and get
(−iω + vik + Dk 2 )C = 0 (6.7)
which in turn leads to
ω = vk − iDk 2 (6.8)
The first term is well-known from e.g. electromagnetic waves (photons) or me-
chanical waves in crystals (phonons), which have linear dispersion relations (for
phonons only for small k). The second term is special for diffusion.
134
where Gn is the Gibbs free energy at position n and F is the external force against
which the motor has to work. Thus equilibrium dictates
k+
= e−β(∆G+F a) (6.11)
k−
using the equilibrium condition from equation 6.11. Thus we get a finite free
velocity at F = 0, then a convex up decay to the stall force Fs and finally a
plateau at negative values (compare 6.1B). Indeed such a force-velocity relation
is known from many motors, e.g. for myosin II (although this is a non-processive
motor, so this is the average result when working in a group) and to some extent
for kinesin.
An alternative scenario would be that the force dependence resides completely in
k− . We then get
v(F ) = a(k+ − k− (F )) = ak+ 1 − eβ(∆G+F a) (6.13)
This force-velocity relation is convex down (compare Fig. 6.1C) and is similar
to the one measured for myosin V, although the divergence to negative values at
large F is of course unrealistic. Fig. 6.1D shows some examples for measured
force-velocity curves and demonstrates that we were able to capture their general
features well with our simple one-state model.
[AT P ]
∆Gh = ∆G0 − kB T ln (6.14)
[ADP ][Pi ]
The first term represents the energetic part of breaking the high-energy bond
in ATP and gives a value around −12.5kB T (to avoid entropic effects, here we
consider very high concentrations, namely M ). The second term represents the
135
(A) (B)
(C) (D)
Figure 6.2: ATP dependence. (A) When the ATP-dependence is only in the forward
rate, then only the free energy barrier height ∆G+ changes when ATP concentration is
changed. (B) When the ATP-dependence is only in the backward rate, then only the
free energy barrier height ∆G− changes when ATP concentration is changed. (C) The
experimental results for kinesin show the linear dependence at low ATP and the plateau
at high ATP predicted by the theory. (D) Force dependence of kinesin for different ATP
concentrations. All four graphs taken from the book Physical Biology of the Cell, chapter
16 on molecular motors.
136
entropic part and corresponds to the law of mass action. For physiological condi-
tions ([AT P ] = mM , [ADP ] = 10µM , [Pi ] = mM ) and the reference concentra-
tion of M to make the argument dimensionless, we get −11.5kB T . Thus together
we have ∆Gh = −24kB T . Note that an ATP-molecule is an energy currency
that is valid twice as much inside the cell than with the reference concentrations,
because the cell keeps ATP at a much higher concentration than ADP. In general,
the free energy gain from ATP-hydrolysis depends on environmental conditions
but usually is between 20kB T and 30kB T . This is usually more than enough for
a molecular motor to perform its powerstroke. With a powerstroke distance of
around 8nm and a stall force of around 5pN (typical values for kinesin), we have
an energy of 40nmpN ≈ 10kB T , which correspond to an efficiency of around 0.5,
if one ATP-molecule gives around 20kB T .
Like for the force dependence, the equilibrium considerations do not completely
determine the ATP-dependence of the jump rates. We again consider the two
extreme cases that the external factor affects only one of the two rates. We first
consider that ATP only affects the forward rate. We now use Kramers reaction
rate theory that states that the transition rate k depends on attempt frequency
Γ and barrier height ∆G as
k = Γe−β∆G (6.15)
The exponential dependence between barrier height and transition rate is also
known as Arrhenius factor in physical chemistry and should not be understood
to be a Boltzmann factor. This law means that the transition rate goes down
dramatically (exponentially) if the barrier height increases.
For our problem we can write
to relate the two barrier heights to each other (we count the two barrier heights
as positive, while the free energy difference is negative, therefore the minus sign
on the left). If we assume that only the forward rate is changed by ATP, then this
means that only ∆G+ is changed when changing ATP (compare Fig. 6.2(A)).
We now can write the two rates as
[AT P ]
v = a(k+ ([AT P ]) − k− ) = a(Γ+ e−β(∆G− +∆G0 ) − Γ− e−β∆G− ) (6.19)
[ADP ][Pi ]
where except for [AT P ], all other quantities are constant. Thus it increases
linearly with ATP-concentration. However, this relation cannot be valid at high
[AT P ], because then the barrier disappears (the left well is pushed up over the
barrier) and Kramers theory is not valid anymore. Thus this must be a result for
low [AT P ].
137
B A
A B
B A
A B
Figure 6.3: (A) In the two-state model, the motor in state 0 has to convert to state 1.
During this process, it can remain stationary or take a step to the left. (B) The motor
in state 1 has to convert to state 0. During this process, it can remain stationary or take
a step to the right.
As the second case, we assume that only the backward rate is ATP-dependent.
Now only the barrier height ∆G− is assumed to be ATP-dependent (compare
Fig. 6.2(B)) and we get
k+ = Γ+ e−β∆G+ (6.20)
−β∆G− −β(∆G+ −∆Gh )
k− = Γ− e = Γ− e (6.21)
and therefore
[ADP ][Pi ]
v = a(k+ − k− ([AT P ])) = a(Γ+ e−β∆G+ − Γ− e−β(∆G+ −∆G0 ) ) (6.22)
[AT P ]
Thus now the dependence is inverse in [AT P ]. The divergence at low [AT P ]
cannot be valid because then the barrier vanishes (the right well is pushed up
over the barrier). Thus this result says that the dependence should plateau at
high [AT P ].
Together, we now have found that the velocity should increase linearly at low
[AT P ] and the plateau at a constant value at high [AT P ]. This is exactly the
experimentally measured dependence for all motors. Fig. 6.2(C) shows this for
kinesin. The plateau velocity is typically around µm/s and the crossover concen-
tration at sub-mM . Fig. 6.2(D) shows the force-velocity relation for kinesin for
different ATP concentrations.
138
move to the left and right only through the state 0 and 1, respectively, with the
rates given in Fig. 6.3. The corresponding master equations are
dp0 (n, t) + − − +
= kA p1 (n − 1, t) + kB p1 (n, t) − kA p0 (n, t) − kB p0 (n, t) ,(6.23)
dt
dp1 (n, t) − + + −
= kA p1 (n + 1, t) + kB p0 (n, t) − kA p1 (n, t) − kB p1 (n, t) .(6.24)
dt
These dynamical equations can be solved by using a continuum limit and a
Fourier ansatz. To obtain the drift velocity, however, it is sufficient to use steady
state arguments. We introduce the total probabilities to be in state 0 or 1:
P
Pi (t) = n pi (n, t). The different positions in the master equation now do not
matter anymore because we sum over them. The dynamic equations for the total
probabilities follow from above as
dP0 (t) + − − +
= (kA + kB )P1 (t) − (kA + kB )P0 (t) , (6.25)
dt
dP1 (t) − + + −
= (kA + kB )P0 (t) − (kA + kB )P1 (t) . (6.26)
dt
These linear equations can be solved easily. For the steady state, however, we do
not even have to do this, but we simply set the time derivatives to zero and get
+ − − +
(kA + kB )P1 (t) = (kA + kB )P0 (t) (6.27)
from which we can also read of the effective rates defined by v = a(k+ − k− ):
+ +
(kA kB )
k+ = − − + + , (6.31)
(kA + kB + kA + kB )
− −
(kA kA )
k− = − − + + . (6.32)
(kA + kB + kA + kB )
We note that this form of the overall rate is rather generic and also follow e.g.
for the steady state approximation of Michaelis-Menten kinetics for enzymatic
processes. The enumerator is a product of rates because it describes a sequence
of steps, and the denominator is a sum of all rates, which is the speed limit for
the process.
139
Figure 6.4: Schematics of the myosin II motor cycle with five steps. The two states
above the line are unbound and the three below are bound. Binding and unbinding
corresponds to establishing and removing a crossbridge between motor and actin filament,
respectively. Myosin II makes two powerstrokes with the second one being a catch bond
(slower under force), such that muscle can work also under high load. It also has a safety
exit (the slip path) to unbind if force is too large. It finally unbinds from the rigor state
by binding ATP. If no ATP is present, the system crosslinks and becomes rigid (this
happens in dead bodies).
140
A B
C D
Figure 6.5: (A) Typical potentials for the two-state isothermal ratchet model. Both the
potentials Wi and the transition rates ωi have to be periodic. Motion ensues if detailed
balance for the rates is broken. (B) Motion corresponds to a tilt in the effective potential
Wef f . (C) The simplest example would be the switch between a flat potential and
an asymmetric sawtooth potential. (D) Two sawtooth potentials with a shift between
them as well as with localized transition rates are more effective. This scheme seems to
correspond to the hand-over-hand motion of kinesin. (A) and (B) from Julicher review,
(C) and (D) from Duke review.
We start with the Fokker-Planck equation from eq. 6.4 and rewrite it as continuity
equation
ṗ(x, t) + J ′ (x, t) = 0 (6.33)
with the flux
J = vp − Dp′ (6.34)
We next assume overdamped dynamics (no mass) and write the velocity as
v = µF = µ −W ′ − Fext (6.35)
where the mobility µ is the inverse of the friction coefficient and the force F is
divided into force from a potential W (x) and an external force Fext (the minus for
Fext means that the force is positive if it acts to the left, against the movement of
the motor, like above for the one-state motor). We also make use of the Stokes-
Einstein relation D = kB T µ, which is an example of the fluctuation-dissipation
theorem (D is fluctuation, µ is dissipation, and the two are not independent of
each other, but related by temperature). Thus we now have for the flux
J = −µ (W ′ + Fext )p + kB T p′ (6.36)
Together Eqs. 6.33 and 6.36 define the FPE as we use it here.
141
We now write the Fokker-Planck equations for two states with switching between
them:
Note that the switching terms with the rates ω1 and ω2 differ only by a minus sign
between the two states. For simplicity, in the following we assume µ1 = µ2 = µ.
The two potentials Wi define the two states. An extreme case would be a sawtooth
potential and a flat potential, for example because one state is charged and the
other is not. At any rates, the potentials Wi (x) and the switching rates ωi (x)
should have the same periodicity with unit cell size l, because this is imposed
by the track with repeat distance l. We also note that the two potentials are of
purely physical (passive) origin and therefore should not be tilted, what means
that
∆Wi = Wi (x = l) − Wi (x = 0) = 0 (6.40)
Otherwise the motor would exhibit motion to the right simply because it moves
down a gradient. What we want to model here is the opposite, namely the
fact that molecular motors spontaneous generate motion in a non-tilted energy
landscape by locally burning energy, without a global gradient.
The beauty of this model is that we do not have to specify potentials and rates
to get the general results we are after. We define total probability and total flux
as X X
P (x) = pi (x), J(x) = Ji (x) (6.41)
i i
If we sum up the first two equations from Eq. 6.39, the switching terms drop out
and we get a FPE without source terms:
We now calculate the total flux. Using pi = λi P and the product rule we get
X
J = −µ λi Wi′ + kB T λ′i + λi Fext )P + kB T λi P ′ (6.44)
i
" #
X
= −µ ( λi Wi′ + 0 + Fext )P + kB T P ′
(6.45)
i
Thus we have exactly the general form for the flux, compare eq. 6.36, if we define
an effective potential by
ˆ x !
X
Wef f (x) = dx′ λi (x′ )Wi′ (x′ ) (6.46)
0 i
142
Note that the switching rates ωi enter indirectly through the occupancies λi .
We now consider the case without external force Fext and ask under which con-
ditions the system can generate directed motion by itself. Such motion appears
if the effective potential is tilted over one period. We therefore define
!
ˆ l X
∆Wef f = dx′ λi (x′ )Wi′ (x′ ) (6.47)
0 i
and ask under which conditions this quantity becomes finite. We then immedi-
ately see that two conditions have to be satisfied:
• The potentials Wi (x) and/or the transition rates ωi (x) have to be asymmet-
ric under x → −x. Otherwise the integrand was symmetric and the integral
vanished. A simple example for this would be an asymmetric sawtooth po-
tential. Then the transition rates in principle could be symmetric, but one
can show that this is not very efficient, so one expects both potentials and
transition rates to be asymmetric.
• The switching rates have to break detailed balance, which means that the
system has to be out of equilibrium. Otherwise the steady state distribution
would be the Boltzmann distribution
e−Wi /kB T
λi (x) = P −W /k T (6.48)
ie
i B
Thus the integrand in Eq. 6.47 would be a total derivative and the integral
would vanish.
These two conclusions are non-trivial and must be valid for any specific motor
model. One also can show that they are true for the case µ1 =
6 µ2 .
For many purposes, it is useful to define the deviation from equilibrium. This
can be done by writing
143
A B F
Figure 6.6: (A) We consider an ensemble of motors that is coupled through a rigid back-
bone. Each motor can bind to the filament and slide down the potential corresponding to
the current state. (B) For sufficient deviation Ω > Ωc from equilibrium, the force-velocity
relation becomes negative and spontaneous symmetry breaking with finite velocities v+
and v− occurs at cero forcing. As the external force is varied, a hysteresis loop emerges.
(C) Spontaneous motion occurs because excitation causes a dip in p1 (x) that then moves
to the right with velocity v. This effectively increases the numbers of motors pulling
further to the right. Like during a phase transition, this is an instability.
144
potentials with periodicity l, compare Fig. 6.6(A). We consider the variable x
to be cyclic, thus we only have to deal with a unit cell with 0 ≤ x ≤ l. We
again consider the Fokker-Planck equation for two states, but now for a common
velocity v:
Normalization reads
1 1
p1 (x, t) + p2 (x, t) = ⇒ p2 = − p1 (6.55)
l l
and dx(p1 + p2 ) = 1.
´
We now consider steady state, ṗi = 0. Together with the normalization, the first
of the two Fokker-Planck equations now gives
ω2
v∂x p1 = −(ω1 + ω2 )p1 + . (6.56)
l
The force balance (or momentum conservation) gives
ˆ l
v v
Fext = − Fint = + dx p1 ∂x (W1 − W2 ) (6.57)
µ µ 0
where the constant term drops out because we integrate over ∂x W2 and W2 is
periodic.
For specific choices of the potentials Wi and the rates ωi , these equations for
the steady state p1 (x) can now be solved. This will then lead to a force-velocity
relation Fext (v). Here we want to proceed with generic properties of this theory
and therefore make a Taylor expansion in small velocity v:
∞
X (n)
p1 (x) = p1 (x)v n . (6.58)
n=0
145
with
ˆ
(n) (n)
F = dx p1 ∂x (W1 − W2 ) . (6.61)
For simplicity we next specify for symmetric potentials (Wi (x) = Wi (−x)), so all
even coefficients vanish (F (0) = F (2) = · · · = 0) and the force-velocity relation
Fext (v) becomes anti-symmetric:
1
Fext = ( + F (1) )v + F (3) v 3 + O(v 5 ) (6.62)
µ
thus this quantity is positive and the only solution to the force-velocity relation
at cero forcing (Fext = 0) is v = 0. Thus with detailed balance, no spontaneous
motion can occur. However, at Ωc > 0 the coefficient F (1) can be negative with
F (1) = −1/µ. Then finite values for v become possible for Ω > Ωc , compare Fig.
6.6(B), and the velocity rises as ±(Ω − Ωc )1/2 . Thus for sufficiently large devi-
ation from equilibrium, the system spontaneously starts to move. The scaling
exponent 1/2 is typical for the mean field theory and the system has to sponta-
neously break symmetry to move either right or left with velocities v+ and v− ,
respectively. Note that in contrast to the single motor case, spontaneous motion
ensues even for symmetric potentials; for single motors, the asymmetric potential
is required to give it its direction, but for multiple motors, the system is persistent
and a spontaneous symmetry break occurs. If one now switches on the external
force, one can move the velocity away from its value at the transition point,
compare Fig. 6.6(B). For example, if the ensembles moves to the right with ve-
locity v+ , one can pull it to smaller velocities with a negative external force Fext .
However, at a critical value of Fext , this branch looses stability and jumps to a
negative velocity. The same works in the other direction and there is a hysteresis
loop. In general, the mathematical structure of this theory is exactly the same
as for the second order phase transition of the Ising model for ferromagnetism.
Velocity v corresponds to the magnetization M , external force Fext to the mag-
netic field H, and the deviation from equilibrium Ω to the inverse temperature
β. If the system works against an external spring, oscillations occur, as observed
often in experiments with molecular motors. A famous example are spontaneous
oscillations of hear bundles in the inner ear, which lead to otoacoustic emissions.
Fig. 6.6(C) shows the main mechanism generating the instability leading to
spontaneous motion. As the motors are excited at the minimum of W1 (x), one
gets a dip in p1 (x). This dip moves to the right, effectively repopulating the
motor population pulling to the right. Thus any fluctuation will be increased
and the system is unstable. The same mechanism is at work in phase transition,
when the system does not counteract a fluctuation.
146
6.8 Master equation approach for motor ensembles
The main advantage of the ratchet models is that they allow us to analyze the
fundamental requirements for motion. In order to describe the function of mo-
tor ensembles in close comparision to experiments, however, one usually turns
to master equations. Here we discuss a simple version that describe coopera-
tive cargo transport by an ensemble of motors2 . Rather than focusing on the
spatial position of the motor ensemble, we rather will later enforce movement
in one direction, but focus on the physical limits of this movement, that is on
the possibility that the walk stops because the ensemble looses contact with its
track. Therefore we now will consider the relevant internal state of the ensemble,
namely the number of bound motors. Similar approaches have been used before
to describe the internal dynamics of adhesion clusters3 .
We consider N motors, of which 0 ≤ n ≤ N are bound at any time t. The
variable n(t) is described by a one-step master equation:
ṗn = ǫn+1 pn+1 + πn−1 pn−1 − (ǫn + πn )pn (6.64)
where ǫn and πn are dissociation and association rates, respectively. The station-
ary state leads to the detailed balance condition
ǫn+1 pn+1 = πn pn (6.65)
and this allows us to calculate the steady state probabilities in a recursive manner:
n−1
Y πi
pn = p 0 (6.66)
i=0
ǫi+1
PN
where normalization n=0 pn = 1 gives us the starting condition
N n−1
!−1 N n
−1 Y
!−1
X Y πi X πi
p0 = 1+ = 1+ . (6.67)
n=1 i=0
ǫi+1 ǫ
n=0 i=0 i+1
From here we define a few quantities of interest. The average number of bound
motors is
N
X pn
Nb = n (6.68)
n=1
1 − p0
where we have excluded the state n = 0 from the sum and have normalized in
respect to the bound states only. The average velocity is
N
X pn
vef f = vn (6.69)
n=1
1 − p0
2
Our treatment is taken from Stefan Klumpp and Reinhard Lipowksy, Coooperative cargo
transport by several molecular motors, PNAS 201: 17284-17289, 2005. Compare also the related
paper by Melanie Müller, Stefan Klumpp and Reinhard Lipowsky, Tug-of-war as a cooperative
mechanism for bidirectional cargo transport by molecular motors, PNAS 105: 4609-4614, 2008,
which generalizes this ansatz to two competing motor ensembles pulling in opposite directions.
3
Compare Thorsten Erdmann and Ulrich S. Schwarz, Stability of adhesion clusters under
constant force, Phys. Rev. Lett., 92:108102, 2004.
147
A
B C
Figure 6.7: (A) Definition of the one-step master equation for cooperative transport by
motor ensembles. (B) Distribution of walking distances for group size N = 1, 2, 3, 4, 5
kinesin motors without load. For large N , these distributions become very flat and their
averages grow. (C) Force-velocity relation for N = 1, 2, 3, 5, 10 kinesin motors. The
effective stall force increases and the curve changes from linear to concave-down. From
Klumpp and Lipowsky PNAS 2005.
ǫef f (1 − p0 ) = π0 p0 . (6.70)
The inverse of this would be the average time < ∆tb > it takes for the ensemble
to unbind. We can calculate
p0 π0 ǫ0
ǫef f = π0 = P = (6.71)
1 − p0 N −1 Qn πi
1+
PN −1 Qn πi
n=0 i=0 ǫi+1 n=1 i=1 ǫi+1
Finally we define the average walking distance < ∆xb >. This can be simply
achieved by replacing ǫn by ǫn /vn and πn by πn /vn in the formula for < ∆tb >
(thus replacing inverse time steps by inverse step sizes):
N n
−1 Y
!
v1 X πi vi+1
< ∆xb >= 1+ (6.72)
ǫ1 ǫ v
n=1 i=1 i+1 i
148
6.8.1 Without load
Our model definition is now concluded and to continue, we have to specify rates
and velocity for the different states n. We first consider the case of vanishing
external load. We set
ǫn = nǫ, πn = (N − n)π, vn = v (6.73)
assuming that each bond dissociates and associates with constant rates ǫ and π,
respectively, independent of the others, thus leading to the combinatorial factors,
and that velocity v is independent of state. We define γ = π/ǫ, the dimensionless
(re)binding rate. The probability distribution now is simply a binomial distribu-
tion, because each bond is open and closed with the probabilities 1/(1 + γ) and
γ/(1 + γ), respectively:
! n N −n !
N 1 γ N γn
pn = = (6.74)
n 1+γ 1+γ n (1 + γ)N
With some work, one can check that this agrees with the general formulae given
above. In particular, we have p0 = 1/(1 + γ)N .
The average number of bound motors follows from the average of the binomial
distribution (with the normalization to the bound states):
1 1 γ γ(1 + γ)N −1 γ
Nb = < n >= N= N≈ N (6.75)
1 − p0 1 − p0 1 + γ (1 + γ) − 1
N 1+γ
with the last expression being valid in the limit of very large N . In the limit of
large γ we get Nb ≈ N . The effective unbinding rate follows as
p0 N γǫ
ǫef f = π0 = (6.76)
1 − p0 (1 + γ)N − 1
and therefore the average bound time is
1 (1 + γ)N − 1
< ∆tb >= = (6.77)
ǫef f N γǫ
From here we get the average run length
v
< ∆xb >= v < ∆tb >= [(1 + γ)N − 1] (6.78)
N γǫ
Thus we see that it increases exponentially with N , thus larger clusters can walk
for much longer distances as long as γ > 1. In the limit of very weak binding, we
make a Taylor expansion in γ and find
v (N − 1)
< ∆xb >≈ [1 + γ] (6.79)
ǫ 2
The first term corresponds to the single motor. In this case, the additional
increase due to cooperativity is only linear in N .
In the case of kinesin, we have v = µm/s, ǫ = 1 Hz and π = 5 Hz, thus γ = 5
and run length increases exponentially with motor number. For N = 5, we are
already up from 1 to 311 µm, and 100 motors give one meter. However, this is
only a statement on the average. One can also calculate the full distribution and
finds that it becomes very broad for large collectives.
149
6.8.2 With load
In order to deal with the case of mechanical load F , we use the linearized force-
velocity relation:
F
vn (F ) = v(1 − ) (6.80)
nFs
Very importantly, here we account for the fact that force F is distributed over the
bound motors (load sharing), thus dissipating its effect over the cluster. While
the association is usually assumed to be force-independent, for the dissociation
we have to take force into account:
The second equation (Bell-relation) takes into account that dissociation is expo-
nentially increased by force, as explained by Kramers theory. Here again we also
take load sharing into account. For kinesin, the stall force Fs and the detachment
force Fd are 6 and 3 pN, respectively. If one now evaluates the formulae given
above for these rates, one finds that increasing N leads to a much slower decay in
the force-velocity relation, and changes its character from linear to concave-down.
150
Chapter 7
Diffusion
∂
ˆ ˆ ˆ ˆ
ρdV = ~
ρ̇dV = − (ρ~v )dA = − ∇(ρ~ ~ v )dV
∂t
~ v ) = 0 ρ=const
⇒ ρ̇ + ∇(ρ~ ⇒ ∇ ~ · ~v = 0 (7.1)
151
Figure 7.1: For many molecules moving together, we can assume a smooth flow field
~v (~r, t).
where f~e is the external force (e.g. gravity) and f~i is the internal force (viscous
force and pressure). For our further discussion we neglect f~e and focus on the
internal forces, that follow from the stress tensor σij = η(∂i vj + ∂j vi ) − pδij with
shear viscosity η as
f~i = ∇
~ · σ = η△~v − ∇p
~
∂~v ~ · ~v △t + · · ·
d~v (~r, t) = ~v (~r + ~v dt, t + dt) − ~v (~r, t) = △t + ~v · ∇
∂t
Under the assumption that the fluid is incompressible, we thus get a non-linear
partial differential equation of motion, the Navier-Stokes equation:
∂~v ~ v = η△~v − ∇p
~
ρ + (~v · ∇)~ Navier-Stokes equation (7.2)
∂t
| {z }
inertial forces
152
Together with the continuity equation, we have four equations for the four vari-
ables ~v (~r, t) and p(~r, t).
In the Navier-Stokes equation inertial and viscous forces compete with a relative
magnitude called Reynolds number:
ρv 2
L ρvL
ηv = =: Re Reynolds number (7.3)
L2
η
g m
m
• For humans (swimming) Re = cm3 s
10−3 P a·s
= 106
g µm
µm
• for molecules/cells Re = cm3 s
10−3 P a·s
= 10−6
The difference in the Reynolds numbers is twelve orders of magnitude, which em-
phasizes the different character of the motion. For humans swimming in water the
predominant forces are of inertial origin, but for the molecules the hydrodynamic
forces have viscous character.
We can say that molecules live in an Aristotelian world, while our intuition comes
form a Newtonian one. Being a molecule or cell means living in syrup — every-
thing is very viscous. Contrary to our every day intuition, the movement of
molecules in water stops if there is no force to sustain it. In this limit (Re ≪ 1)
the Navier-Stokes equation becomes the Stokes equation
~
η△~v = ∇p (7.4)
For example, from the Stokes equation one can calculate the Stokes force required
to drag a sphere of radius R with velocity v through a viscous fluid (η), see figure
7.2:
FS = 6πηRv (7.5)
We call the constant term 6πηR the friction coefficient ξ.
The most drastic difference between the Navier-Stokes and the Stokes equations
is that the time-dependence has dropped out. This implies that if we stop and
reverse a motion, we get back exactly to the initial state. Therefore a microscopic
swimmer cannot move by reversing a one-dimensional process, like the up and
down movement of a rigid appendix (scallop theorem, because the macroscopic
scallop moves by opening and closing its shells). Microswimmers have evolved
two main strategies to beat the scallop theorem: either a flexible appendix (like
the eukaryotic flagellum of e.g. sperm cells or green algae, making the forward
and backward strokes very different) or a helical appendix that can turn without
reversing its motion (like the bacterial flagellum).
Despite the viscous nature of its environment, there is one type of forces, which
always keeps molecules on the move, namely thermal forces. The resulting motion
153
Figure 7.2: A sphere with radius R moving through water under the action of force FS .
with
∆x2
D= diffusions coefficient
2∆t
Note that the mean squared displacement grows only linearly with time, like
the mean squared end-to-end distance of polymers grows only linearly with their
contour length.
We next want to relate the diffusion constant to the viscosity. We start with
some simple arguments and later present a rigorous derivation. We now consider
a small biomolecule of mass 100kDa:
1 kB T
mv 2 = equipartition
2 2
1 1
kB T 2 4.1 pN · nm 2 m
⇒v = = = 10
m 100 kDa s
This is the typical velocity between collisions. Dynamically molecular velocity
can be thought to be generated by some thermal force F , acting between the
collisions:
mẍ = F
F
⇒ ∆x = (∆t)2
2m
∆x F ∆t F
⇒v = = =
∆t 2m ξ
154
Here ξ is the friction coefficient.
2m
ξ =
∆t
∆x2
⇒ Dξ = m 2 = mv 2 = kB T
∆t
kB T Einstein relation
⇒ D= connecting diffusion and friction (7.7)
ξ
kB T
D= Stokes-Einstein relation (7.8)
6πηR
Thus the more viscous the environment, the smaller the diffusion constant.
A numerical example for a biomolecule is
4.1 pN nm −6 cm
2 (10 µm)2
D= = 10 = (7.9)
6π10−3 Pa · s 1 nm s s
2 · 10−6 cms
2
2D
∆x = = = 0.2 Å
v 10 ms
and
∆x 0.2 Å
∆t = = = 2 ps
v 10 ms
One can conclude that mean free path length and collision time are extremely
small.
We now present a more rigorous derivation of the Einstein relation starting from
the Fokker-Planck (or Smoluchowski) equation:
→
− → − −
ṗ = − ∇ · J (→
r)
~ ~
J(~r) = −D∇p + p · ~v (7.10)
| {z } |{z}
diffusion drift
F~ = −∇U
~ = ξ~v
155
→
−
and look how the system reacts. In equilibrium ṗ = 0, which means that ∇ ·
→
− →
J (−
r ) = 0. If the system is confined, the flux has to vanish. That means
~
∇U
~
D∇p = p · ~v = −p
ξ
U U
⇒ p(~x) = p0 exp − = p0 exp −
Dξ kB T
| {z }
Boltzmann factor
kB T kB T
⇒ D= = Stokes-Einstein relation
ξ 6πηR
We again get the Stokes-Einstein relation for the diffusion coefficient. In contrast
to the scaling argument above, this derivation is rigorous.
By this method intense laser light is used to bleach (destroy) the fluorophores in
a well-defined region, e.g. a stripe with diameter 2a on a plate with length 2L,
see figure 7.3.
For example one can monitor the concentration of fluorophores over time at a
given point in space, e.g. c(x = 0, t). Then the recovery kinetics allows one to fit
the diffusion constant D.
For the full calculation we have to know the initial conditions:
156
Figure 7.3: FRAP. Concentration of fluoreophores on a plate of length 2L. a) System
before bleaching. b) System at time t = 0, when the laser light bleaches a stripe of width
2a. c) Bleached spot d) Recovery of the system after some time ∆t.
c0 c0
c0
increasing t
c0
c 2 c 2
0 0
-L -a 0 a L -L -a 0 a L
x x
Figure 7.4: Fluorescence intensity as a function of the space coordinate after bleaching.
On the left side you can see the bleached spot at the initial moment. On the right-hand
side you can follow the recovery of the intensity with time.
157
c0
−L < x < −a
c(x, 0) = 0 −a < x < a
c a<x<L
0
We also need to define some boundary conditions, e.g. no flux at the boundaries:
c′ (x, t) |x=±L = 0
Now we have an even function with the right boundary conditions, and we still
have to find the coefficients A0 and An .
The cosine functions form an orthogonal system:
ˆL
xnπ xmπ
cos cos dx = Lδnm
L L
−L
∞ ∞
!
n2 π 2
∂A0 X ∂An xnπ X xnπ
⇒ + cos = D −An 2 cos
∂t n=1
∂t L n=1
L L
∂A0
⇒ = 0
∂t
∂An Dn2 π 2
= − An
∂t L2 !
n2 π 2
An (t) = An (0) exp −(D 2 )t
L
We still need the initial values (Fourier transform of the step function):
ˆL
1 L−a
A0 (0) = dx c(x, 0) = c0 = c∞
2L L
−L
ˆL
1 xnπ
An (0) = dx c(x, 0) cos
L L
−L
Here we can use the symmetry of the system and calculate the integral just for
x between a and L and multiply the result by 2.
158
c0
c0
2
c
0
t
ˆL
2c0 xnπ
= dx cos
L L
a
2c0 L xnπ L
= sin a
L nπ L
2c0 anπ
= − sin
nπ L
Thus we get
" ∞
! #
n2 π 2
a 2 X 1 anπ xnπ
c(x, t) = c0 1− − sin exp −(D 2 )t cos
L π n=1 n L L L
This sum converges slowly, which means that if one wants a precise result, one has
2
to take many modes into account. Each mode has its own time scale tn = DπL2 n2 .
This makes the calculations very difficult. A convenient choice for simplification
of the equation is a = L2 . And to qualitatively understand the result we look just
at the first mode, see figure 7.5, e.g.
x = 0
n = 1
1 2 t
c(x = 0, t) = c0 − exp −
2 π t1
159
∆n
t
Figure 7.6: Photon count over time n(t) for molecules freely diffusing in and out of the
focal volume.
function 2 :
τD τD τD
F (t) = exp − J0 + J1 (7.11)
2t 2t 2t
2
Here F (t) is normalized average over the bleached spot; τD = wD , with w the
radius of the focused beam and J0 and J1 are the Bessel functions of zeroth and
first order.
The FCS method was developed and published in the 1970s3 . N fluorescent
molecules pass through the sub-femtoliter detection volume V of a confocal mi-
croscope (the standard instrument for FCS). One can measure the number of
photons n(t) in the photomultiplier as a function of time, see figure 7.6.
Then one calculates the auto-correlation function
hδn(t)δn(t + τ )i
G(τ ) =
< n >2
where δn(t) = n(t)− < n > are the variations of n. The auto-correlation function
is the convolution of the signal at time t with the signal at the same place after
2
Axelrod, D., Koppel, D. E., Schlessinger, J., Elson, E. & Webb, W. W. Mobility measure-
ment by analysis of fluorescence photobleaching recovery kinetics. Biophys. J. 16, 1055-1069
(1976); Soumpasis, D. M. 1983. Theoretical analysis of fluorescence photo- bleaching recov-
ery experiments. Biophys. J. 41:95-97; Sprague, B.L. et al., Analysis of binding reactions by
fluorescence recovery after photobleaching, Biophys. J. 86: 3473-3495 (2004).
3
D. Magde, E. Elson, and W. W. Webb, Phys. Rev. Lett., 29,705 (1972); D. Magde, E. L.
Elson, and W. W. Webb, Biopolymers, 13,29 (1974); S. R. Aragon and R. Pecora, J. Chem.
Phys. 64, 179 (1976)
160
n or m alized var ian ce
GHΤL
G2
Figure 7.7: Auto-correlation function. This fit measures the self-similarities of the signal
after a lag time τ .
z2
!
2(x2 + y 2 + K2
)
ρ(x, y, z) = ρ0 exp −
w2
w is the e−2 size of the detection volume. For the laser power we get
3/2
π
ˆ
P = d~rρ(~r) = ρ0 w3 K
2
| {z }
effective detection volume V
which means that < n >=< c > V = N is the number of molecules in the focus.
ˆ
⇒ δn(t) = d~rρ(~r)δc(~r, t)
161
with the boundary condition δc(~r = ∞, t) = 0.
Now, we transform to Fourier space:
∂(δc(~k, t))
= −Dk 2 δc(~k, t)
∂t
⇒ δc(~k, t) = δc(~k, 0) exp −Dk 2 t ,
We first evaluate
ˆ
< δc(~r, 0)δc(r~′ , τ ) >= d~k exp i~k r~′ < δc(~r, 0) δc(~k, τ ) >
| {z }
=δc(~k,0) exp(−Dk2 τ )
1
ˆ ˆ
< δc(~r, 0)δc(r~′ , τ ) >= d~k exp(i~k r~′ ) exp(−Dk 2 τ ) ~ exp(−i~k r”)
dr” ~ < δc(~r, 0)δc(r”,
~ 0) >
(2π)3 | {z }
=<c>δ(~ ~
r −r”)
Thus we assume that spatial correlation only exists at the same position and that
the second moment is proportional to the average (typical for a Poisson process).
Therefore
<c>
ˆ
< δc(~r, 0)δc(r~′ , τ ) >= d~k exp(−i~k(~r − r~′ )) exp(−Dk 2 τ )
(2π)3
Therefore
3 !
(~r − r~′ )2
<c> π 2
< δc(~r, 0)δc(r~′ , τ ) > = exp −
(2π)3 Dτ 4Dτ
The integral separates in the three spatial dimensions. We calculate it for the
162
x-direction:
! ! !
2x2 2x′2 (x − x′ )2
ˆ ˆ
′
dx dx exp − 2 exp − 2 exp −
w w 4Dτ
! ! !
(u2 + 2ux′ + x′2 ) x′2 −u2
ˆ ˆ
u=x−x′ ′
= du dx exp −2 exp −2 exp
w2 w2 4Dτ
ˆ !
2 1 4x′2 4ux′
ˆ
2 ′
= du exp − + u dx exp − exp −
w2 4Dτ w2 w2
| {z }
1
πw2 2 u2
= exp
4 w2
!1 ˆ
πw2
2
1 1 2
= du exp − 2
+ u
4 w 4Dτ
!1 !1
πw2 2
π 2
= 1 1
4 w2
+ 4Dτ
!1
πw 1 2
1 w2
= (4Dτ ) 2 with τD =
2 1 + ττD 4D
3 2 !1
< c > ρ0 π 2 πw 1 3 πKw 1 2
⇒ G(τ ) = (4Dτ ) 2
τ
3
(2π) 2 N 2 Dτ 2 1+ τ 2 1+ K 2 τD
τD
! !1
<c>V =N ρ20 1 1 2
=
N 1 + ττD 1+ τ
KτD
1 1 w2
G(τ ) = with τD = (7.14)
N τ
3 4D
1+ τD
2
Thus a fit to the experimental data would allow us not only to extract D, but
also N (w is known). This calculation has been extended in many ways, e.g.
reactions, 2D rotational diffusion, 2D diffusion (membranes), anomalous diffusion
of polymers, protein complexes, etc.
163
Figure 7.8: Two examples of diffusion to capture processes. On the left-hand side is a
scheme of actin filament growth. The actin monomers diffuse in the cytosol and when
they meet the growing filament, they attach to it. On the right hand side the binding
of a transcription factor to the DNA double helix is depicted. At first the transcription
factor binds unspecifically and slides along the DNA molecule, until it finds the specific
binding side.
receptors is 104 . The ligands undergo Brownian diffusion and obey the diffusion
equation
ċ = D△c (7.15)
For a well-mixed situation we would have standard reaction kinetics
kf
R+L⇋C
kr
where kf is the forward reaction constant and kr the reverse reaction constant.
ċ = kf RL − kr C
ċ=0 RL kr
⇒ = = KD
C kf
Here KD is the dissociation constant as specified by the law of mass action. For
low affinity receptors a typical value for KD is 10−6 M , and for high affinity ones
KD = 10−12 M .
Note that R, C are not concentrations, but absolute numbers in this case. To
understand the KD we look at the half-occupancy
L = KD
⇒R = C
164
Figure 7.9: Cell with receptors R on its surface. The cell explores its immediate envi-
ronment, by responding to the signal from those receptors, when ligands L are binding
or unbinding. Before binding to the receptors the signal molecules first have to diffuse
to the surface of the cell.
where RT is the total number of receptors, we see that this is again a Poisson
process. If R = 104 and L = KD , then the concentration fluctuations δc c = 1%,
which is a relatively large number. One percent deviation of c can lead to cellular
response and change in the cell behaviour (e.g. chemotaxis by E. coli). The cell
has to integrate the signal in order to react adequately to the changes in the
environment. To understand this we have to investigate the role of the diffusion.
In order to understand the effect of diffusion, we consider the cell to be a perfect
spherical adsorber of radius R. The diffusion equation in spherical coordinates
reads
1
ċ = D 2 ∂r (r2 ∂r c)
r
In steady state (ċ = 0) ,r2 ∂r c must be constant. For the perfect spherical adsorber
c(R) = 0
R
c(r) = c0 1 −
r
For r = R the concentration is zero and for r → ∞ it acquires some constant
value c0 , see figure 7.10.
The local concentration flux is j = −D∂r c. Then the total current on the sphere
is
2 R 2
J = 4πR j(R) = −4πR D c0 2 = −4πRDc0
R
165
c0
R
c
R
r
Figure 7.10: Relaxation of the ligand concentration c(r) over time for a perfect adsorbing
sphere. At the surface of the sphere the concentration is zero, and far away it reaches a
constant value c0 .
Thus J depends linearly on R because the gradient and the surface area give
counteracting trends. c0 is the driving force here, because adding more ligands
results in stronger current.
The diffusion determines the rate at which the ligands hit the surface of the cell.
The association rate for pure diffusion is
|J|
k+ = = 4πDR
c0
This famous result is called the Smoluchowski rate. Using typical values for pro-
teins (R = 4nm, D = 2(10µm)2 /s), one gets a rate of k+ = 6×109 1/(M s), which
is 3-4 order of magnitude larger than experimentally measured diffusion-limited
association rates4 . The reason is the assumption of isotropic reactivity. Solc
and Stockmayer have shown that one gets the experimentally measured values if
one takes anisotropic reaction patches and rotational diffusion into account. We
conclude that the localization of binding to surface patches effectively decreases
association by 3-4 orders of magnitude. Association can be increased in the range
of the Smoluchowski rate by attractive interactions (steering) between the reac-
tion partners, with the highest known rates around 1010 1/(M s). In figure 7.11
we give an overview over the resulting scenario.
In order to calculate the dissociation rate, we consider a sphere with concentration
c0 on its surface and a perfect sink at infinity:
c0 R
c(r) =
r
R
⇒J = 4πDR2 Dc0 = 4πDRc0
R2
which is basically the same current as for association, but in the reverse direction,
4
For an excellent review on this subject, compare Gideon Schreiber, Gilad Haran, and Huan-
Xiang Zhou, Fundamental aspects of protein-protein association kinetics, Chemical Reviews
109:839-860, 2009
166
Solc & Stockmayer Smoluchowski
enhanced by
electrosta6c
a7rac6on
(reac6on)
167
c0
R
c
R
r
Figure 7.12: Sphere with concentration c0 on the surface. Now the concentration of
ligands decays from its initial value to zero away from the sphere.
168
Diffusion Electrostatics
concentration ∆c = 0 potential ∆φ = 0
~
flux ~j = −D∇c electrical field E~ = −∇φ
~
~ 1 ~ A~
total current J = ~jdA charge Q =
´ ´
4π Ed
Figure 7.13: Model receptors as small patches with radius s on the sphere with radial
dimension R. Now the flux is not rotationally symmetric, because the current flows only
towards the patches. That is why there is no exact solution in this case. But if we look
at a slightly bigger radius, R + d, then the flux is again uniform.
From the relations given in table 7.1 one finds that J = 4πDRc0 as calculated
above. Using the same analogies, one can find e.g. that the flux from a halfspace
onto a disc of radius R is J = 4DRc0 (same scaling, but no factor π).
Until now we have treated the cell as uniformly reactive. In practice, it is covered
only by N receptors, which can be modeled as N disc-like adsorbers of capture
radius s ≪ R. The rest of the sphere is considered to be reflecting. For simplicity
we consider the diffusion-limited case. We use the electrostatic analogy, where
the total current J is equal to the ratio of driving force (voltage) over the Ohm’s
resistance R0 :
c0
J=
R0
Because in this case there is no radial symmetry, see figure 7.13, there is no
longer an exact solution. For a good approximation (5% compared to computer
simulations) we introduce a slightly larger radius R + d and look at r > R + d,
where the flux is uniform again. Since d ≪ R by definition, we put d → 0 in the
calculations. For better illustration we introduce the equivalent electric current
diagram, see figure 7.14 , where R0R is the Ohm’s resistance of the sphere and
169
Figure 7.14: Electrical current schema equivalent to the model of receptor patches on the
cell surface. It represents a combination of N receptors with resistance R0s connected in
parallel.
k +c 0
N12
N
Figure 7.15: Current flow as a function of the number of receptors N . For small N the
dependency is linear, and with higher N the current J gets to saturation.
1
R0R =
4πDR
1
R0s =
4Ds
R0s
and R0 = R0R +
N
1
⇒J = c0 1 1
4πDR + 4DsN
(
4πDRc0 4πDRc0 for N → ∞
= = (7.16)
1 + πR
sN 4DsN c0 for N → 0
170
cal number for surface receptors. The adsorbing fraction of the sphere is then
πs2
only N
4πR2
= 1.6 · 10−4 and the distance between neighbouring receptors is only
1
4πR2
N = 140 nm = 140 · s. This shows that diffusion is not efficient for trav-
2
eling long distances, but is very efficient for exploring space on short distances.
Because random walks are locally exhaustive, few receptors are sufficient. This
explains why the cell can accommodate so many different systems in parallel:
they do not need to be space filling and leave sufficient space for other systems.
171
172
Chapter 8
Reaction kinetics
173
e-‐,O2
H20
H+
ΔV
ADP
H+ ATP
Figure 8.1: The electron transport chain in the inner membrane of mitochondria ends
with cytochrome C oxidase (left), that builds up a proton gradient (∆V = 200 mV). This
gradient is then used by the ATP synthase (right) to convert ADP into ATP, the energy
currency of the cell. Typical concentrations for ATP and ADP in the cell are mM and
10 µM, respectively.
where ∆V is the membrane potential around 200 meV storing the en-
ergy (compare figure 8.1). This gradient is then used by ATP synthase
to produce ATP molecules. The network motifs of metabolic networks are
strongly shaped by optimality in regard to the use of chemical and energy
resources.
174
TF
RNAp transcription translation
mRNA Protein
ribosome
DNA
Gene
Ø degradation Ø
Figure 8.2: A gene on the DNA is expressed by the sequential processes of transcription
and translation. The resulting protein can feedback on this or another gene through
transcription factors. Feedback onto itself can be either positive (gene is turned on)
or negative (homeostasis). Feedback onto other genes can be either upregulating or
inhibiting. In this way, a complex gene expression network emerges. Humans have around
24.000 genes resulting in more than 60.000 proteins (through alternative splicing).
175
GF
GFR
MAPKKK MAPKKK-‐P
MAPKK MAPKK-‐P
MAPK MAPK-‐P
TF
Figure 8.3: One of the most important signal transduction pathways is the MAPK-
pathway (mitogen activated protein kinases, leads to cell division). A growth factor
activates a GF-receptor, e.g. epidermal growth factor (EGF). Then a cascade of three
levels of kinases is used to amplify the signal, eventually leading to activation of a tran-
scription factor (TF). The hierarchical structure of the pathway leads to ultrasensitivity.
In general, many signaling pathways are branched and interact with other signaling path-
ways. In this way, a signal transduction network emerges. Again such networks based
on protein-protein interactions are stored in databases. The interaction strengths can be
measured e.g. with the yeast-two-hybrid method.
176
8.2 Law of mass action
The mathematics to describe biochemical networks starts with the law of mass
action. To understand this law we consider a container filled with reactive species
of molecules with homogenous concentration (well-mixed assumption, good for
small reaction volumes such as an E. Coli cell). We want to know how fast the
reactions take place and how fast the concentrations of the reacting molecules
change. The answer is given by the law of mass action. This law states that
the reaction speed is proportional to the number of combinations of reacting
molecules. The number of combinations is in turn proportional to the concen-
trations of the reactants. The proportionality constant is called reaction rate
constant k, which summarizes all effects going in the effective reaction (move-
ment in space, encounter, final complexation).
Consider for example two species of substrates (S1 , S2 ), which react with each
other to produce two identical product molecules (P ). Schematically the reaction
looks like this
k+
S1 + S2 ⇋ 2P
k−
Here k− and k+ are the rate constants for the backward and forward reactions,
they are called dissociation and association rates. The forward reaction depends
on two processes. First the molecule have to diffuse through the medium and get
close to each other and then they can react. These two situations are described
by the rate constants kd (diffusion) and ka (association). These we can combine,
to get the rate constant for the forward reaction:
1 1 1
= + (8.3)
k+ kd ka
For the forward reaction the number of combinations of S1 and S2 is NS1 ,S2 =
S1 S2 ∼ s1 s2 , where s1 , s2 are the concentrations of the two substrates. The speed
of the forward reaction is
v+ = k+ s1 s2
The number of combinations of P for the reverse reaction is
f or P ≫1
NP = P (P − 1) ≈ P 2 ∼ p2
where p is the concentration of the product molecules. The speed of the reverse
reaction is
v− = k− p2
The net speed of the reaction can be calculated as the sum of the speeds of the
forward and the reverse reactions, taken with the appropriate sign:
v = v+ − v− = k+ s1 s2 − k− p2 (8.4)
177
describe the rates at which those concentrations change. The rate at which the
concentration of substrate S1 changes with time is given by
ds1 ds2
= −v+ + v− = −v = −k− s1 s2 + k− p2 =
dt dt
The rate, at which the concentration of S2 changes in time, is the same as the
rate for S1 . For the change of the concentration of the product over time we get
dp
= 2v+ − 2v− = 2v = 2(k+ s1 s2 − k− p2 )
dt
We should keep in mind that the total number of molecules is conserved
s1 + s2 + p = const
This means that we can get rid of one variable and work with one ODE less.
1
Note that v has the dimension M s and k+ and k− have the dimension M ·s .
Now we want to look at the set of ODE for the concentrations for stationary
states. We assume a system that does not change in time, so that ds ds2
dt = dt = 0
1
and dp
dt = 0. For the speed this means v = 0, so the forward and the reverse
reactions occur with the same rate v+ = v− . Keq is the rate constant that
describes the equilibrium state between forward and backward reactions and is
called the equilibrium constant:
k+ (peq )2
Keq := = eq eq (8.5)
k− s1 s2
This is the law of mass action for this special case.
For a generalized description we have to look at M substrate species with mul-
tiplicity µi , and N product species with multiplicity νi . The reaction can be
written as
M
X k+
N
X
µi Si ⇋ νi Pi
k−
i=1 i=1
For the reaction speed we get a generalized expression
M
Y µi
v+ = k+ Si
i=1
YN
v− = k− Piνi
i=1
Q
k+ N
Piνi
⇒ Keq = = Qi=1
M µi (8.6)
k− i=1 Si
178
In this section we used a classical reaction kinetics approach, which assumes
deterministic macroscopic systems with large number of molecules, spatial ho-
mogeneity and fast diffusion. This approach cannot describe systems with small
number of molecules, with slow movements (e.g. in a crowded environment) or
with spatial heterogeneities. For these situations, more sophisticated approaches
exist, in particular particle-based computer simulations.
k2
EX1 + X ⇋ EX2
k−2
Q is called binding polynomial, it is the sum over all possible ligation states. Q
can be interpreted as a partition sum.
We are interested in the fraction of all E molecules in the solution, which have i
bonded ligands. We look at the cases for i = 0, 1, 2:
1
po =
Q
k1eq x
p1 =
Q
eq eq 2
k1 k2 x
p2 =
Q
179
Figure 8.4: Many enzymes can bind to more than one substrate. If the binding is
cooperative, than the affinity of the remaining binding sites changes. The fraction of
molecules with no bonded substrates is depicted in blue, the fraction with one bonded
molecule is represented by the red line and the fraction with two bonded molecules is the
green line, all represented as functions of the substrate concentration x.
We can see in figure 8.4 how the different states are populated in respect to the
ligand concentration x.
We can also estimate the number of bound molecules of type X
k1 x + k1 k2 x2
< i >=
Q
To illustrate the effect of cooperativity we put in some typical values and look at
the distribution, assuming p0 = p2 :
k1 = k2 = 100
⇒ xmid = 0.01 and Q = 3
1
⇒ p0 = p 1 = p 2 =
Q
k1 = 103 ≫ k2 = 10
⇒ xmid = 0.01 and Q = 12
1 10
⇒ p0 = p 2 = and p1 =
Q Q
In such cases the intermediate states are more probable.
• For activating binding
k1 = 1 ≪ k2 = 104
⇒ xmid = 0.01 and Q = 2.01
1 0.01
⇒ p 0 = p2 = and p1 =
Q Q
180
This indicates that strong reinforcement of binding depopulates the intermediate
state.
For kk21 → 0
n RLn
R + nL → C, KD = = const (8.13)
C
where KD has been defined somehow differently in order to get the same dimen-
sion as L. We now get
C Ln
pb = = n (8.14)
R+C KD + Ln
thus n is the Hill coefficient of this system and KD the threshold value.
In order to obtain a more mechanistic understanding for n = 1, we now consider
a lattice model and treat it with the canonical formalism. Consider a lattice with
Ω sites. We have L ≪ Ω ligands which are distributed over the lattice, each with
unbound energy ǫu . The corresponding partition sum is
Ω! ΩL −βLǫu
Zu = e−βLǫu ≈ e (8.15)
L!(Ω − L)! L!
181
We now add one receptor to the system that can bind a ligand with binding
energy ǫb . This gives a second part to the partition sum:
Zb (L/Ω)e−β∆ǫ
pb = = (8.17)
Zu + Zb 1 + (L/Ω)e−β∆ǫ
182
and Bela Novak, Sniffers, buzzers, toggles and blinkers: dynamics of regulatory
and signaling pathways in the cell, Current Opinion in Cell Biology 15:221-231,
2003. The corresponding plots are taken from this paper, too. Very similar
conclusions about the high level functions of network motifs have been reached
also in engineering (control theory).
183
Synthesis
and
degrada.on
Phosphoryla.on
(buzzer)
sigmoidal
Adapta.on
(sniffer)
Figure 8.5: Three simple networks motifs that all three lead to stable fixed points. Left:
the reaction scheme. Middle: plot of gain and loss terms. The intersection is the fixed
point. Because loss / gain dominates at higher / lower values, it is stable. Right:
the steady state or bifurcation diagram as a function of input signal S. The hyperbolic
response of the buzzer becomes sigmoidal if cooperativity or enzyme kinetics are involved.
For the sniffer, we show the time course, because only in this way the adaptation response
becomes clear (the steady state for R does not change as a function of S, but it shows a
transient each time S is changed).
184
Posi%ve
feedback
by
mutual
ac%va%on
Figure 8.6: Two examples of positive feedback both leading to bistability and toggle
switches. An essential element is the non-linearity of the system leading to the S-curves.
To get more complex responses, we can now add a second species. We discuss
first a feed forward motif that leads to adaptation:
Ṙ = k1 S − k2 XR, Ẋ = k3 S − k4 X (8.28)
We now come to the first example with feedback. A typical example for positive
feedback would be mutual activation, such that the protein R upregulates a
185
Nega%ve
feedback
without
delay
(homeostasis)
Nega%ve
feedback
with
delay
(oscilla%ons)
Figure 8.7: Negative feedback can lead either to homeostasis or oscillations, depending
on whether the system has sustained delay or not.
Ṙ = k0 + k1 S − k2 R − k3 E(R) (8.31)
186
to achieve a constant power output. Here we study this for a biochemical system
(chemostat).
We consider a species R that inhibits a species E, which however helps to produce
R. Thus if too much / too little of R exists, less / more is produced (production
on demand). The dynamical equation is
Ṙ = k0 E(R) − k1 SR (8.32)
For steam engines, it was found that the delay between the sensor and the actu-
ator can lead to oscillations. For engines, this has to be avoided, but in biology,
oscillations are often put to good use, e.g. to anticipate the break of dawn (cir-
cadian rhythm). In order to get some delay into the negative feedback cycle, we
introduce a new species. We consider X upregulating Yp and Yp upregulating Rp ,
but Rp inhibiting X. Thus we have a negative feedback as before, but it takes
more time to propagate the signal. The following system of equations describe
such a system:
Ẋ = k0 + k1 S − k2 X − k7 Rp X (8.33)
k3 X(YT − Yp ) k4 Yp
Ẏp = − (8.34)
Km3 + (YT − Yp ) Km4 + Yp
k5 Yp (RT − Rp ) k6 Rp
Ṙp = − (8.35)
Km5 + (RT − Rp ) Km6 + Rp
(8.36)
For all three species as a function of S, there exists a region in which oscillations
occur. The amplitude of these oscillations appears and vanishes in a smooth
manner (Hopf-bifurcation). Although a steady state does not exist in this re-
gion, the system still can be represented in a bifurcation diagram by showing
the minimal and maximal values. Biochemical oscillations often use the negative
feedback motif (blinker). Often they also involve space (waves resulting from
reaction-diffusion systems) and / or mechanics (like the otoacoustic oscillations
in our ears). Hopf bifurcations are very common in biological systems which want
to regulate the amplitude of the oscillations. If the oscillation would jump to a
finite amplitude, this could lead to a catastrophe (this can occur e.g. for bridges
suddenly collapsing after a small change in the bifurcation parameter).
187
Figure 8.8: Energy landscape of a reaction with (dark red line) and without (light red
line) enzymes. Enzymes do not take part of the reaction itself, they rather decrease the
energy barrier ∆E of the system, thus enabeling the substrates to easily convert into
product molecules.
k1 k
S + E ⇋ SE −→
2
P +E (8.37)
k−1
188
Here k1 is the rate constant of formation of the enzyme-substrate complex, k−1 is
the rate constant of dissociation of the enzyme-substrate complex, and k2 is the
catalysis rate constant. In effect the substrates are turned into products. The
classical analysis of this central scheme of enzyme kinetics goes back to Leonor
Michaelis und Maud Menten in 1913. A more rigorous analysis has been worked
out later by different mathematicians, most prominently by Lee Segel. For more
details, see chapter 6 of the book by Murray on mathematical biology.
From now on we will use the following abbreviations for the concentrations of the
reacting molecules: s = [S], e = [E], c = [SE], p = [P ]. To define the system of
the molecules, we need to define the initial conditions for the concentrations
We get the ODE describing the concentration changes of the molecules in time
by applying the law of mass action on our reaction:
ds
= −k1 es + k−1 c
dt
de
= −k1 es + (k−1 + k2 )c
dt
dc
= k1 es − (k−1 + k2 )c
dt
dp
= k2 c (8.38)
dt
Since the ODE for the product rate is decoupled from the other three, we can
directly give the concentration dependency on time
ˆt
p(t) = k2 c(t′ )dt′ (8.39)
0
ds
=−k1 e0 s+(k1 s+k−1 )c with s(0) = s0
dt
dc
=k
(8.40)
dt 1 e0 s−(k1 s+k−1 +k2 )c with c(0) = 0
We first note that the reaction is irreversible and that all substrate will be used
up eventually. The classical Michaelis-Menten analysis assumes that enzyme is
present in a small quantity and quickly converted into the complex, which then
189
is in a steady state until all substrate has been used up; then also the complex
has to vanish again. Therefore their main assumption is
dc
≈0 (8.41)
dt
The ODE for c now becomes an algebraic relation
e0 s(t)
c(t) = (8.42)
s(t) + Km
and the ODE for s has the form
ds vmax s(t) dp
=− =− (8.43)
dt s(t) + Km dt
where we have defined the Michaelis constant
k−1 + k2
Km := Michaelis-constant
k1
and the maximal velocity
vmax = k2 e0
The famous result for the production rate therefore reads
dp vmax s
vp = dt = s+Km (8.44)
thus c(t) is not constant as assumed and in particular, it does not satisfy the
initial condition c(0) = 0. This shows that this analysis is not consistent and
needs to improved. In particular, we have to understand in which sense it might
be relevant after all.
190
Figure 8.9: Michaelis-Menten kinetics. The reaction rate as a function of the substrate
concentration, growing linearly with s, for small s, and than asymptotically approaching
vmax . Source wikipedia.org
191
Figure 8.10: This diagram depicts the behaviour of the substrate u(τ ) and SE complex
v(τ ) as a function of the nondimensional time. One can also see the progression of the
concentration of the unbound enzyme ee0 . Source: Mathematical biology book from J.D
Murray
du0
= −u0 + (u0 + K − λ)v0
dτ
0 = u0 − (u0 + K)v0
u0
⇒ v0 =
u0 + K
du0 λu0
= −
dτ u0 + K
⇒ u0 (τ ) + K ln[u0 (τ )] = 1 − λτ (8.50)
Thus this procedure gives us exactly the same result like the standard treatment.
Although this solution could be improved systematically by going to higher order,
this does not solve our problem that somehow our treatment becomes inconsistent
again.
192
The solution to our problem is that we realize that in the perturbation scheme
we make the strong assumption that the solutions are analytical for all times.
However, this is not true for small times. We solve this problem by going from
regular to singular perturbation theory. We now magnify the scale near the zero
by introducing a new time scale τǫ = σ. This allows us to look closely at the
region around τ = 0, because even small changes in τ have large effects on the σ
scale. We use the substitutions u(τ, ǫ) = U (σ, ǫ) and v(τ, ǫ) = V (σ, ǫ) to get
dU
= −ǫU + ǫ(U + K − λ)V
dσ
dV
= U − (U + K)V
dσ
We again approximate U (σ, ǫ) and V (σ, ǫ) by a Taylor expansion in ǫ :
∞
X
U (σ, ǫ) = ǫn Un (σ)
n=0
X∞
V (σ, ǫ) = ǫn Vn (σ) (8.51)
n=0
193
introduction into non-linear dynamics, we refer to the book by Steven Strogatz,
"Nonlinear Dynamics and Chaos", Westview Press, December 2000.
In general, dynamical systems are described by a dynamical (differential) equation
of the state vector ~x = (x1 , . . . , xn ):
f1 (~x)
˙~x = f~(~x) = .. (8.52)
.
fn (~x)
where the overdots denote differentiation with respect to the time t. The vari-
ables x1 , . . . , xn might represent any species whose dynamics we are interested in,
such as concentrations or population levels. The functions f1 (~x), . . . , f2 (~x) are
determined by the specific problem which we want to analyze.
The most general results of NLD are the following:
• For n = 2, oscillations become possible. There are also a few more types of
possible bifurcations.
mẍ + kx = 0
!
0 1
⇒ ~x˙ = ~x = A~x (8.53)
−ω 2 0
Thus, we end up with a linear version of our general dynamic equation. We also
conclude that the harmonic oscillator is in fact a two-dimensional system. This
is an example that oscillations become possible only in two dimensions.
The solutions of equation 8.52 can be visualized as trajectories flowing through
an n-dimensional phase space with coordinates x1 , . . . , xn . In general, non-linear
dynamics deals with analyzing these flow properties, often in a graphical manner.
194
Figure 8.11: Phase portrait of an arbitrary function ẋ = f (x). The pink arrows denote
the direction of flow on a line. Fixed points (ẋ = 0) in one dimension can either be stable
(solid dot) or unstable (hollow dot), meaning that already small perturbations will result
in a flow away from the unstable fixed point.
Also note the intrinsic time scale of the system, τ = 1/f ′(x∗ ).
Examples:
1. Logistic growth
The most prominent model of mathematical biology is the logistic growth
model of a population size N (t). This was suggested to describe the growth
of human populations by Verhulst in 1838:
N
Ṅ = rN (1 − ) (8.56)
K
where N ≥ 0. Here, r denotes the growth rate of the population and K is
the carrying capacity which accounts for population-limiting effects such as
195
(a) (b)
Figure 8.12: a.) Phase portrait of the logistic growth model (equation 8.56). The carrying
capacity K denotes a stable fixed point. b.) The population always approaches K for
large times. Note the qualitatively different behaviors for different initial values of N .
2rN ∗
f ′(N ∗ ) = r − = −r > 0 X
K
Hence, the carrying capacity K is approached by the population from any
initial value N > 0, but in a qualitatively different manner (compare figure
8.12b, calculation not shown here).
2. Autocatalysis
An autocatalysis is a chemical reaction which exhibits a positive feedback
to produce the species X:
k1
A + X ⇋ 2X (8.57)
k2
With the law of mass action we can immediately identify the differential
equation for the concentration of X:
ẋ = k1 ax − k−1 x2 (8.58)
Thus we find that this equation is essentially the same as for logistic growth.
For a flow on a line, the trajectories can only increase, decrease or stay constant
(compare figure 8.11), but they cannot cross themselves. Hence, no oscillations
in a 1D system are possible. There are other ways to arrive at this conclusion.
We can understand the first-order differential equation as an overdamped one-
dimensional system and define a potential V (x) such that
dV
ẋ = f (x) = − (8.59)
dx
196
(a) ẋ = −x (b) V (x) = x2 (c) qualitative time depen-
dence
Figure 8.13: Phase portrait, potential and time evolution of a linear f (x) (upper panel)
and a cubic f (x) (lower panel).
This integral always exists. We can now compute the derivative of V with respect
to the time t
dV 2
dV dV dx
= · =− ≤0 (8.60)
dt dx dt dx
which implies that the system loses energy until it has reached a local minimum,
i.e. a stable fixed point. To illustrate this, two examples are shown in figure
8.13. At any case, because the system comes to a halt, it cannot oscillate. This
is different for n = 2, where energy is lost in one form but can be converted into
another form.
8.7.2 Bifurcations
At first glance, the flow on the line seems to be a rather boring problem: With
no oscillations, trajectories always end up in a stable fixed point or head out
to infinity. However, as a function of model parameters, fixed points can be
created or destroyed, or their stability can change. These qualitative changes in
the system dynamics are called bifurcations and the parameter values at which
they occur are so-called bifurcation points. Mathematical bifurcations usually
represent some transition, instability or threshold in the physical model.
Bifurcations can be visualized in a bifurcation diagram which plots the fixed
points x∗i against the parameter r which is varied. Stable fixed points are usually
depicted as solid curves whereas unstable fixed points are depicted as dotted lines.
Note, that there can be several parameters ri , leading to a multidimensional
bifurcation diagram which then becomes increasingly complicated.
For n = 1, there are four types of bifurcations.
197
(a) ẋ = r + x2 with r < 0 (b) ẋ = r + x2 with r > 0 (c) bifurcation diagram
The saddle-node bifurcation is the basic mechanism by which fixed points are
created or destroyed. A typical example is shown in figure 8.14. With in-
creasing r, the parabola is shifted in the positive ẋ-direction and the two fixed
points x∗1 and x∗2 move closer to each other until they annihilate when r = 0.
For r > 0, there is no fixed point at all and therefore no stable solution of the
dynamic equation. Hence, the bifurcation point is located at r = 0 (compare the
bifurcation diagram in figure 8.14c).
198
(a) ẋ = rx − x3 with r < 0 (b) ẋ = rx − x3 with r > 0 (c) bifurcation diagram
199
Figure 8.18: Phase plane trajectory of the harmonic oscillator. The closed orbit cor-
responds to a periodic motion. The amplitude of the periodic motion depends on the
initial condition.
ẋ v
=
v̇ −ω 2 x
⇒ −ω 2 x dx = v dv
1 2
⇒ 2v + 12 ω 2 x2 = const (8.62)
1 1 2
⇒ mv 2 + kx = const (8.63)
2
| {z } |2 {z }
kinetic energy potential energy
We thus find that the system follows an elliptic path in the phase plane (equation
8.62, compare figure 8.18). Obviously the system is not overdamped as it is
necessarily the case in 1D, since the energy of the system is conserved for all
times (equation 8.63).
In biology, one often encounters oscillations. One out of many examples is the
circadian rhythm which can be found in every single mammalian cell. However,
such systems typically are non-linear in contrast to the linear harmonic oscillator,
because linear systems have some disadvantages:
200
2. Perturbations stay small, but they are not corrected.
The flow-pattern we have observed for the harmonic oscillator is only one out
of several flow patterns occurring in n = 2. Linear systems are the reference
point for non-linear ones. We will now introduce a systematic classification of
flow patterns in linear systems, and then expand this idea to non-linear ones:
!
a b
ẋ = A · ~x with A = (8.64)
c d
where λ1/2 are the eigenvalues of A and ~v1/2 the corresponding eigenvectors. The
eigenvalues are determined with the characteristic polynomial:
(a − λ)(d − λ) − cb = 0
⇒ λ2 − (a + d)λ + (ad − cb) = 0
| {z } | {z }
=:τ =trA =:∆=det A
√
τ ± τ 2 −4∆
⇒ λ1/2 = 2 (8.66)
For 4∆ > τ , the eigenvalues are complex and the system oscillates. Hence,
∆ = τ 2 /4 defines the boundary of the region where oscillations occur. Further
analysis will result in six possible flow patterns which are illustrated in a phase
diagram (figure 8.19).
As already mentioned, we can use our results of the linear system as a reference
point for non-linear systems in a linear stability analysis similar to the one-
dimensional case. Recall the dynamical equation:
! !
ẋ f (x, y)
=
ẏ g(x, y)
Well-behaved non-linear systems can be expanded around their fixed points (x∗ , y ∗ ):
one first investigates the the two nullclines defined by ẋ = f (x, y) = 0 and
ẏ = g(x, y) = 0 (figure 8.20). Secondly, one investigates their intersections which
are the fixed points. We perform a linear stability analysis around them:
u := x − x∗ , v = y − y∗
201
Figure 8.19: Different flow patterns for a linear 2D system. The gray area marks the
area where oscillations occur in the system. Note that the center is unstable in regard
to small perturbations.
⇒ u̇ = ẋ = f (x, y) = f (x∗ + u, y ∗ + v)
= f (x∗ , y ∗ ) + ∂x f (x∗ , y ∗ )u + ∂y f (x∗ , y ∗ )v + O(u2 , uv, v 2 )
| {z }
0
! ! !
similar for v u̇ ∂x f ∂y f u
⇒ = · (8.67)
v̇ ∂x g ∂y g v
| {z }
=:A=J(x∗ ,y ∗ )
Now one can analyze the matrix A according to our procedure of the linear system
above and thus classify the fixed points. Usually this is a good starting point for
a full analysis of the non-linear system.
202
(a) (b)
(c)
Figure 8.21: Van der Pol oscillator. (a) Flow in the (x, y)-phase plane: the system slowly
follows the cubic nullcline and then quickly zips to the other side. (b) Flow in the (x, v)-
phase plane. (c) Approach to the stable limit cycle from outside and inside. In both
cases, we eventually end up with the sawtooth oscillations typical for the van der Pol
oscillator.
Different from the harmonic oscillator, this one has negative damping at small
amplitude and positive damping at large amplitude. Thus the oscillator is driven
up by energy input if it relaxes and it is dampened if it is agitated. From this,
an intermediate state of sustained oscillators is emerging with an amplitude that
is independent of initial conditions.
We analyze this system in the regime µ ≫ 1 according to Lienard. We rewrite it
in the following way:
d
− x = ẍ + µ(x2 − 1)ẋ = (ẋ + µF (x)) := ẇ (8.69)
dt
where we defined F (x) = (x3 /3 − x). We therefore have a new set of dynamical
equations
ẋ = w − µF (x), ẇ = −x (8.70)
We now define y := w/µ and therefore get
203
(a) (b)
Figure 8.22: a.) A positive feedback loop in signal transduction. The signal S will be our
bifurcation parameter. The response protein R acts as a kinase for the inactive enzyme E,
i.e. it adds an inorganic phosphate to the enzyme to produce the active enzyme Ep . The
latter can enhance the production of R. b.) Signal-response curve for the steady state
response R∗ in the form of a one-dimensional bifurcation diagram exhibiting bistability.
In this example, bistability forms a reversible switch.
These equations are easy to understand: the nullclines are y = F (x) and x = 0
and the steady state at (0, 0) is not stable because the Jacobian reads
!
−µ(x2 − 1) µ
A= (8.72)
−1/µ 0
thus we have τ = µ and ∆ = 1 at (0, 0). For small µ, we have an unstable spiral,
and for large µ, we have an unstable node.
The flow in the (x, y)-plane can be understood as follows (compare figure 8.21).
We start with an initial condition y−F (x) = O(1). Then the velocity is fast in the
x-direction (O(µ)) but slow in the y-direction (O(1/µ)). Thus it zips to the right
side until it hits the cubic nullcline. If y − F (x) = O(1/µ2 ), then both velocities
become comparably small (O(1/µ)). The system now crosses the nullcline and
slides along it until it reaches its extremum. Then it accelerates again in the
horizontal direction and the second half of the cycle begins. Altogether we have
a sequence of four parts (fast, slow, fast, slow). We are always driven to the same
stable limit cycle independent of initial conditions. The resulting sawtooth-like
pattern is typical for a relaxation oscillator, which slowly builds up some tension
and then quickly relaxes it, compare figure 8.21.
204
(a) (b)
Figure 8.23: a.) Open and closed bonds of an adhesion cluster. b.) Possible trajectory
of the number of closed bonds N in time for F = 0.
with2
2uK Goldbeter-Koshland
G(u, v, J, K) = p (8.76)
B+ B2
− 4(v − u)uK function
B = (v − u) + vJ + uK
Though highly non-linear, the main equation (equation 8.73) is now one-dimensional
and its steady-state solution shown in figure 8.22b. We see a window of bistabil-
ity separated by two saddle-node bifurcations. This switch between high and low
responses is a very important control element in biological systems.
with (Nt − N ) denoting the number of open bonds and a force-dependent disso-
F
ciation constant kof f = ko e F0 ·N . Hence, the force couples the bonds and makes
the system non-linear.
2
A Goldbeter and D E Koshland: An amplified sensitivity arising from covalent modification
in biological systems, PNAS 1981
205
(a) (b)
Figure 8.24: a.) Graphical solution of equation 8.80 for large and small f . Note, that
for large f there is no stable solution. b.) Bifurcation diagram of the system. At the
bifurcation point (fc , Nc ) the two fixed points collapse. This is a typical example for a
saddle-node bifurcation.
kon
Ṅ = 0 ⇒ k0 N ∗ = kon (Nt − N ∗ ) ⇒ N ∗ = Nt (8.79)
kon + k0
f
dN
dτ = −e N N + γ(Nt − N )
SS f
⇒ N e N = γ(Nt − N ) (8.80)
fc fc /Nc
γ= e (8.83)
Nt
206
Figure 8.25: Sketch of a positive feedback loop implemented into a gene regulatory
network. The gene is transcribed into mRNA (Y ) which is then translated into the
protein X. The protein can form clusters of size 2 which then act as a transcription
factor enhancing their own gene expression. Both mRNA and protein are subject to
degradation processes within the cell.
fc − Nc
1− =
Nc Nt − Nc
Nt Nc
⇒ fc =
Nt − Nc
fc fc
⇒ = +1 (8.84)
Nc Nt
207
(a) (b)
Figure 8.26: a.) Nullclines of the genetic switch model. The fixed points were numbered
according to the text. The vector field is vertical on the line ẋ = 0 and horizontal on the
line ẏ = 0. The directions can be determined by noting the signs of ẋ and ẏ. b.) Phase
portrait of the genetic switch.
immediately be formulated.
ẋ = y − ax (8.86)
x2
ẏ = − by (8.87)
1 + x2
The nullclines are then given by
y = ax (8.88)
x2
y = (8.89)
b(1 + x2 )
and are shown in figure 8.26a denoting that there are at most three fixed points.
Since the parameter a governs the slope of the curve ẋ = 0, the two fixed points
2 and
3 collide and disappear as a increases. This behavior reminds us strongly
of that of a saddle-node bifurcation in one-dimensional systems.
Now, the positions of the three fixed points are determined:
x∗2
ax∗ =
b(1 + x∗2 )
⇒ ∗
x1 = 0 (8.90)
√
1 ± 1 − 4a2 b2
⇒ x∗2/3 = (8.91)
2ab
The fixed points
3 collide if 1 = 4a2 b2 which defines the critical parameter
2 and
value ac = 1/(2b). At the bifurcation, the fixed point is given by x∗2/3 = 1. Hence,
x∗2 < 1 and x∗3 > 1.
In order to the types of the fixed points, we perform a linear stability analysis
with a Jacobian:
!
−a 1
A = 2x
(1+x2 )2
−b
208
1 we find for the determinant ∆ = ab > 0 and for the trace τ = a + b :
For
,
2x∗
∆ = ab −
(1 + x∗2 )2
x∗2 2ab
with abx∗ = : = ab −
1+x ∗2 1 + x∗2
x∗2 − 1
⇒ ∆ = ab · ∗2 (8.92)
x +1
2 as a saddle point. On the other
Since x∗2 < 0, we get ∆ < 0 which classifies
hand,
3 is a stable node, because with x3 > 1 we now that 0 < ∆ < ab and
∗
hence τ 2 − 4∆ > (a − b)2 > 0. The phase portrait is plotted in figure 8.26b.
If a < aC = 1/(2b), or in other words, mRNA and protein degrade sufficiently
slowly with ab < 1/2, then the system is bistable and acts like a genetic switch.
Depending on the initial conditions, the gene is either "on" or "off". The stable
manifold of the saddle acts like a threshold ("separative"), separating the two
sinks of attraction.
In general, all four types of bifurcation discussed for n = 1 can occur for n = 2.
They keep their one-dimensional nature, with the flow in the extra dimension
being essentially simple attraction or repulsion (example for "center manifold
theory"). The example of the genetic switch shows this for a saddle-node bifur-
cation in n = 2.
Every living cell needs energy in the form of ATP. In oxidation, this is produced
from food sources such as glucose. One molecule of glucose (C6 H12 O6 ) can be
used to produce up to 28 molecules of ATP. The biochemical process of breaking
down sugar to obtain energy is called glycolysis. Since the production of ATP
should adapt to its need, some kind of feedback has to be implemented into the
pathway.
Indeed, glycolytic oscillations have been discovered in 1964 in yeast and muscle
extracts (monitored by fluorescence of NADH, periods are 5 and 20 minutes,
respectively). These oscillations are produced by the phosphofructokinase (PFK),
whose activity is controlled by ATP/ADP (figure 8.27a). A simple mathematical
model goes back to Sel’kov in 19683 . Inhibition by ATP and activation by ADP
makes sure that the pathway becomes active only when needed. However, because
ATP is also a substrate for PFK, this leads to an instability (figure 8.27b).
3
E Sel’kov: "On the Mechanism of Single-Frequency Self-Oscillations in Glycolysis. I. A
Simple Kinetic Model," Eur. J. Biochem. 4(1), 79-86, 1968.
209
(a)
(b)
Figure 8.27: a.) Simplified scheme of the glycolysis pathway in a cell. Note, that ATP
acts both as an inhibitor and an activator for its own production. b.) Simple model for
x and y representing the conflicting role of ATP.
ẋ = −x + ay + x2 y (8.93)
2
ẏ = b − ay − x y (8.94)
at fixed
⇒ ∆ = a + b2 >0 (8.97)
point
b4 + (2a − 1)b2 + (a + a2 ) !
τ =− >0 (8.98)
a + b2
√
⇒ b = 21 (1 − 2a ± 1 − 8a) (8.99)
210
(a) (b) (c)
Figure 8.28: a.) Nullclines of the dynamical equations in the glycolysis oscillator. The
direction of flow was determined by the signs of ẋ and ẏ. b.) Phase diagram of the system.
The gray area marks the values (a, b) which lead to oscillations in the system (unstable
fixed point), whereas outside this area the fixed point is stable and no oscillations occur.
Varying b with a fixed will turn on and off the oscillations (green arrow) which is typical
for a Hopf bifurcation. c.) Sketch of a stable limit cycle. For large t, every flow pattern
in the phase plane will end up in this limit cycle, irrespective of its initial conditions.
Figure 8.29: Sketch of a neuron with intra- and extracellular ionic composition (A−
stands for organic anions like nucleic acids). Signals of other nerve cells are received at
the dendrites (passive cables) and are then passed along the axon (active cable) in form
of an action potential with speed v ≈ 10 − 100 m/s. At the synapses, voltage-gated ion
channels lead to Ca-influx, which in turn leads to vesicle fusion and neutrotransmitter
release. Neurotransmitters act as signals for the post-synaptic nerve cells. If an axon
fires or not is determined at the axon hillcock, where the input signals are summed up.
closed trajectory, a so-called limit cycle (figure 8.28c). Perturbations are coun-
terbalanced and initial conditions do not matter. This motif is often encountered
in biology. Note that in contrast to the bistable switch produced by autocatal-
ysis, here autocatalysis leads to oscillations because there is an extra variable x
that can be used up. Thus the same motif can lead to very different consequences
depending on its context.
211
intracellular (mM) extracellular (mM) Nernst potential (mV)
K+ 155 4 -98
N a+ 12 145 67
Cl− 4 120 -90
Ca2+ 10−4 1.5 130
Table 8.1: Ion concentration and the Nernst potential in a typical mammalian muscle
cell. Since the Nernst potentials of the different ion species differ strongly, this ionic
distribution is an out-of-equilibrium situation. Resting potentials of excitable cells are
in the range of −50 to −90 mV .
1
p= exp − keV
BT
∼c
Z
c1
2)
⇒ = exp − e(Vk1B−V
T
c2
kB T c1 Nernst
⇒ ∆V = ln (8.100)
e } c2
| {z potential
25 mV
The non-equilibrium situation described above cannot be true for a plain lipid
bilayer. In fact, ions have to be transported across the membrane to keep the
system in non-equilibrium. Transport across the membrane can either be passive
4
For a comprehensive overview, compare James Keener and James Sneyd, Mathematical
Physiology, Springer, second edition 2009.
212
along the electrochemical gradient of the ions through open ion channels or active
via so-called pumps. In the plasma membrane, there are numerous such mem-
brane proteins fulfilling these functions. Of particular interest for the formation
of an action potential are voltage-gated ion channels (for K + and N a+ ) and the
Na/K-pump (figure 8.30).
Ion channels can have two conformations: open and closed. The transition be-
tween these two conformations is called gating. For instance, voltage-gated N a+
channels in neurons are closed when the membrane is at the resting potential, but
they open as the membrane potential becomes more positive. The probability of
being open of this two-state system can be described by
e−β·∆ǫ
popen = (8.101)
1 + e−β·∆ǫ
1
∆ǫ = ∆ǫ0 − qV : = ∗ −V ) (8.102)
1+e βq(V
213
(a) (b)
Figure 8.31: a.) Probability to be open for a voltage-gated ion channel following equation
8.102. The threshold potential is denoted at V ∗ . b.) Action potential of a nerve cell. 1:
Threshold for opening of N a+ channels. 2: N a2+ channels open after 0.5 ms and close
2 ms later, leading to a "depolarization" of the membrane. 3: K + channels open after
2 ms and close several ms later, achieving a repolarisation. 4: Overshoot due to open
K + channels ("hyperpolarization"). 5: Resting state. Now the system needs some time
to establish the initial state ("latency"). Similar depolarization waves, but with Ca2+ ,
are used to have all sarcomeres contracted simultaneously in a muscle.
theoretical work on action potentials in the squid giant axon, using the experi-
mental data acquired in the squid summer season of 1949 at Plymouth, when sci-
entific work had started again in post-war England. The Hodgkin-Huxley model
is the most important mathematical model in biology and is still valid today. It
is actually a four-dimensional non-linear dynamics model which they evaluated
on hand-cranked computers. Without knowing about ion channels and pumps,
they postulated exactly the corresponding mathematical structure. Experimen-
tally, they used the new space clamp techniques (figure 8.32a) which allows for
studying an action potential as a function of time without having to consider
spatial effects. A few important historical milestones concerning excitability are
summarized in table 8.2.
We now turn to the electrical properties of the axon, i.e. we express the biological
system by an electric circuit (circuit diagram shown in figure 8.32b):
Recall that ∆Q = C∆U and C/A = µF · cm−2 for a lipid bilayer. The current
components are given by Ii = gi (∆V − Vi ). In the following back of the enve-
lope calculation we will see that the Nernst potentials can indeed assumed to be
214
1952 Famous Hodgkin and Huxley papers: the dynamics of so-called gates
produce temporal changes in conductivity
1960 Richard FitzHugh and later Nagumo et al. independently analyzed a
reduced HH-model with phase plane analysis, leading to the standard
NLD-model for action potentials
1963 Nobel prize for Hodgkin and Huxley (together with John Eccles, who
worked on motorneuron synapses)
1991 Nobel prize for Erwin Neher and Bert Sakmann for the patch clamp
technique: the molecular basis of an action potential could be
demonstrated directly for the first time
2003 Nobel prize for Roderick MacKinnon for his work (Science 1998) on
the structure of the K + channel, which in particular explained why
N a+ ions cannot pass
May
Andrew Huxley dies at the age of 94; after his work on the action
2012
potential, he revolutionized muscle research (the sliding filament
hypothesis from 1954 and the Huxley model for contraction from 1957
could have earned him a second Nobel prize)
(a) (b)
Figure 8.32: a.) Space clamp of the squid giant axon (here: R = 25 µm) where an ar-
bitrary potential can be clamped to the membrane. b.) Circuit diagram of an excitable
membrane. gi denotes the conductance and Vi denotes the Nernst potential of the re-
spective ion species. Hodgkin and Huxley also introduced a third current component
which they called L ("leakage").
215
constant. The charge transfer during depolarization is given by
C
∆Q = ∆V · · (2πRL)
A
µF
= 100 mV · · 2π · 25 µm · L
cm2
e
= 1010 ·L
cm
whereas the total charge inside the axon is
Q = e · c · (πR2 L)
= e · 100 mM · π(25 µm)2 · L
e
= 1015 ·L
cm
The relative change in charge ∆Q/Q is in the order of 10−5 which is a very small
perturbation. Therefore, the Nernst potentials can well assumed to be constant.
In their model, Hodgkin and Huxley also considered some leakage current compo-
nent L. With Ohm’s law (I = Q̇ = C · ∆V˙ = g∆V ) the first dynamical equation
is given by
˙ 1
∆V = − [gN a (∆V − VN a ) + gK (∆V − VK ) + gL (∆V − VL )] (8.103)
C
Note that this simple equation gives a stable fixed point and no interesting dy-
namics is to be expected except if the gi have some dynamics by themselves.
Hodgkin and Huxley realized for the first time that gN a and gK must have very
different dynamics and came up with equations that gave a very good fit to their
data.
At the core of their model is the assumption that some two-state process must
lead to the breakdown of the conductance, suggesting a hyperbolic response. On
the other hand, however, in their voltage clamp data they observed a sigmoidal
response. They concluded that several two-state processes must be combined.
They introduced three dynamic gates m, h and n that together determine the
conductancies:
gN a = 120 m3 h (8.104)
4
gK = 36 n (8.105)
gL = 0.3 = const. (8.106)
Here the exponents have been chosen simply from fits to the data. Each of the
three gates is a two state process, where the transition rate from "closed" to
"open" is given by αi and the transition rate from "open" to "closed" is given by
βi :
ṅ = αn (1 − n) − βn n (8.107)
ṁ = αm (1 − m) − βm m (8.108)
ḣ = αh (1 − h) − βh h (8.109)
216
The αi and βi are phenomenologically determined functions of v = ∆V but have
the functional form of equation 8.102, which today we understand is appropriate
for ion channels:
25 − v
αm = h i
25−v
10 · exp 10 −1
−v
βm = 4 exp
18
7 −v
αh = exp
100 20
1
βh =
30−v
exp 10 +1
10 − v
αn = h i
10−v
100 · exp 10 −1
125 −v
βn = exp
1000 80
With equations 8.103 to 8.109, the HH-model is a NLD-model with the four
variables ∆V , n, m and h. A numerical solution for n, m and h is shown in figure
8.33. The solution for ∆V corresponds to a typical action potential (figure 8.31b)
as observed experimentally. The HH-model was not only a perfect explanation
of all the data they had measured, but it has stayed valid until today.
Before the gate h has recovered to its steady state, an additional stimulus does
not evoke any substantial response, although the potential itself is close to its
resting value. Therefore latency arises. This is another essential observation
explained by the HH-model. Interestingly, this property also explains why the
injection of a steady current can lead to oscillations. Such a current raises the
resting potential above the threshold of an action potential. The axon fires but
then has to recover before it fires again. Another way to make the HH-model
oscillate is to immerse the axon in a bath of high extracellular K + -ions. This
raises the Nernst potential for K + and therefore the resting potential. If it
becomes superthreshold, oscillations start. Later we will see that these kinds of
oscillations are very similar to the ones we obtain from the van der Pol oscillator.
In the 1960s, Richard FitzHugh analyzed several reduced version of the the HH-
model with phase-plane techniques. He started from the observation that the
HH-model has two fast variables, m and v (the Na-channels activate quickly and
the membrane potential changes quickly) and two slow variables n and h (K-
channels are activated slowly and Na-channels are inactivated slowly). He first
analyzed the fast phase plane and found that it shows bistability, which means
that the system can switch from a resting to an excited state. On the long time
scale, n starts to increase (K-channels open) and h starts to decrease (Na-channels
217
(a) (b)
(c)
Figure 8.33: (a) Time dependence of the gates n, m and h during an action potential.
(b) Resulting time dependence of the conductancies. (c) Resulting action potential. ∆V
and m are fast variables, whereas n and h are slow variables. Taken from Keener and
Sneyd.
218
(a) (b)
Figure 8.34: (a) Fast-slow phase plane of the Hodgkin-Huxley model. The FitzHugh-
Nagumo model essentially has the same structure. (b) The resulting action potential.
Taken from Keener and Sneyd.
close). This leads to the disappearance of the excited state through a saddle-node
bifurcation and the system returns to the resting state.
FitzHugh next analyzed the fast-slow phase plane with v and n, which are called
the excitation and recovery variables, respectively. He found that the v-nullcline
has a cubic shape, that there is one fixed point and that the action potential is
emerging as a long excursion away and back to the fixed point guided by the
cubic nullcline, compare figure 8.34. This reduction of the HH-equations then
motivated him to define an even more abstract model that later became to be
known as the FitzHugh-Nagumo model (Nagumo and coworkers built this
model as an electronic circuit and published in 1964). This model assumes two
variables, one slow (w) and one fast (v). The fast (excitation) variable has a
cubic nullcline and the slow (recovery) variable has a linear one. There is a single
intersection which is assumed to be at the origin without loss of generality. Thus
the model equations are
dv
ǫ = v(1 − v)(v − α) − w + Iapp (8.110)
dt
dw
= v − γw (8.111)
dt
where Iapp allows for an externally applied current, ǫ ≪ 1 and 0 < α < 1.
Typical values are ǫ = 0.01, α = 0.1 and γ = 0.5. The phase plane analysis then
show that an excitation to a small value of v leads to a large excursion (action
potential) leading to the steady state (0, 0). If one injects a current Iapp = 0.5,
the fixed point becomes unstable and a stable limit cycle emerges through a
Hopf bifurction, thus the system becomes oscillatory (essentially it becomes a
van der Pol oscillator). Thus this simple model reproduces the main features of
the HH-model.
Interestingly, the HH- and FN-models have many more interesting features if
studied for a time-dependent current Iapp (t). For example, one finds that the
system does not start to spike if the current is increased slowly rather than in
219
Figure 8.35: Extended electrical circuit diagram of an excitable membrane such as the
squid giant axon.
a step-wise fashion. Thus it does not have a fixed threshold but goes super-
threshold only if the current change is fast. Another interesting observation is
that an inhibitory step (negative step function) triggers a spike at the end rather
than at the start of the stimulation period. Thus the direction in which the
current is changed matters.
where r is the differential longitudenal resistance. We next note that the leakage
current decreases the current along the line (Kirchhoff’s balance of currents):
with IT (x, t) having contributions both from the capacitance and the resistance:
dV (x, t)
IT (x, t) = C + gV (8.114)
dt
Together we therefore get in the continuum limit
d2 V (x, t) dV (x, t)
2
= r(C + gV ) (8.115)
dx dt
p
or with the decay length λ = 1/(rg) and the time constant τ = C/g:
d2 V (x, t) dV (x, t)
λ2 2
−τ =V (8.116)
dx dt
220
This is the famous cable equation, a diffusion equation with damping, with a
similar mathematical structures as the Schrödinger or Fokker-Planck equations.
For a constant g, we deal with the passive cable equation that has been derived
in 1855 by Lord Kelvin for the transatlantic cable. This equation applies to the
dendrites.
For the axons, we have to use the active cable equation, where the damping is
replaced by the dynamic equation of the HH-model
d2 V (x, t) dV (x, t) gN a gK gL
λ2 2
−τ = (V − VN a ) + (V − VK ) + (V − VL ) (8.117)
dx dt g g g
and supplemented by the three equations for the conductivities. It has been found
already by Hodgkin and Huxley that such an equation allows for propagating
wave. To mathematically analyze this equation, we can reduce this system to the
bistable cable equation (in dimensionless form):
dV (x, t) d2 V (x, t)
= + f (V ), f (V ) = −V (V − α)(V − 1) (8.118)
dt dx2
which has stable fixed points at V = 0 (resting state) and V = 1 (excited state)
(for 0 < α < 1). One can show that this equation has a traveling front solution
with a tanh-shape. To get a traveling wave like the action potential, one has to
add an relaxation process which brings the excited state back to the resting state.
If one specifies the cable equation for a cylindrical geometry with radius R, one
notes that g and C both scale with R because they increase with the surface area,
and r scales like 1/R2 because it decreases
p with √
the volume. The propagation
velocity therefore scales as λ/τ ∼ g/(rC 2 ) ∼ R, as found experimentally.
This explains why the squid has evolved a giant axon.
The active cable equation can also be extended to two and three dimensions.
Then one not only gets propagating fronts and waves, but also e.g. spirals,
which are self-sustrained and do not propagate a signal. In the context of heart
and muscle biology, such patterns are the sign of a pathophysiological situations
(stroke, seizure).
After having studied the elementary unit of the brain (the neuron and the action
potential) in large detail, we now ask in which way many neurons work together
in the brain, which has around 1011 neurons and 1014 synapses. Because such
a large system cannot be simulated with all its details, we have to find abstract
representations5 . Such abstract models are increasingly simulated in neuromor-
phic hardware in order to emulate the amazing computational capabilities of the
human brain (compare group of Karlheinz Meier at KIP).
5
For a comprehensive introduction, see the books by Wulfram Gerstner and coworkers, Neu-
ronal Dynamics, Cambridge University Press 2014, and by Peter Dayan and L.F. Abbott, The-
oretical Neurosciences, MIT Press 2001.
221
Figure 8.36: (a) Serial blockface scanning electron microscopy image of a part of a rabbit
retina with 5123 voxels (Denk group). (b) Image reconstruction based on a graphical
model with supervoxels (Hamprecht group). From Bjoern Andres, Medical Image Anal-
ysis 16: 796-805 (2012). Such reconstructions should eventually provide a complete map
of the synaptic connections between the neurons (connectome).
We first note that the way a neuron is built, it is clear how one has to connect it up
into a network: dendrites collect the synaptic currents, they are summed up at the
axon hillcock, the axon might fire or not, and the resulting spike it distributed
to the post-synaptic cells. The exact topology of the network depends on the
part of the brain one is interested in and can be taken from experiments imaging
the brain (e.g. with two-photon microscopy, serial blockface scanning electron
microscopy or super-resolution microscopy, compare figure 8.36). We next note
that the shape of an action potential is stereotyped (as described by the HH-
model) and thus what seems to matter most is the timing of the spikes. A simple
model which takes these considerations into account is the leaky integrate-and-fire
model (LIF). For each neuron, we solve a dynamic equation for the membrane
potential:
dV
C = −g(V − E) + I (8.119)
dt
but we disregard spatial aspects and do not include a relaxation variable. The
first term is the ionic (leak) current and the second the input (synaptic) current.
If a synaptic current arises, the voltage starts to rise. In the LIF, a spike is
emitted if the potential crosses a given threshold from below. The potential then
is reset (this replaces the dynamics of the relaxation variable) and can rise again
(after a latency time). In this way, a spike train is generated in response to a
constant input current. Typical values are 0.2 kHz for the frequency and 10 nA
for the current.
As we have discussed above, the HH-model does not imply a fixed threshold, thus
222
the LIF is too simple in this regard. Moreover it results in a constant frequency
under constant stimulation, while experimentally, the system adapts. In order
to take these two issues into account, the adaptive exponential integrate-and-fire
model (AdEx) has been suggested:
dV
C = −g(V − E) + g∆T e(V −ET )/∆T + I (8.120)
dt
dw
τw = a(V − E) + bτw ρ − w (8.121)
dt
The second term pushes the potential up once the value ET is reached. In fact
once this term dominates, the potential diverges in finite time and a spike is
created. After spiking, the potential is reset and latency is implemented, as in
LIF. w is an adaptation variable: it grows with the potential and with the spike
train ρ, but relaxes on the time scale τw . For constant stimulation, the spiking
frequency therefore decreases (adaptation). The time constant τw is of the order
of hundreds of milliseconds. For a = b = ∆T = 0, we essentially recover the
LIF-model.
A complete model for a neural network also needs a synaptic model. In the
current-based synaptic interaction model (CUBA), a spike creates a typical cur-
rent (e.g. the difference between two exponentials) that then is fed into LIF
or AdEx models. Although CUBA-models are a natural choice for point neu-
ron models such as LIF and AdEx, there is also physical justification for a
conductance-based synaptic interaction model (COBA).
Synapses are also the place for learning and memory (plasticity) and therefore
the synapse between neurons i and j gets a weight wij that has its own dynamics.
For example, Hebb’s rule states that synapses increase their weight if they are
used frequently. With these rules, the neural networks is completely defined and
can be trained for a certain function, e.g. automatic detection of handwriting or
license plates.
223
dimensional system of the type:
Ȧ = F (A, B) + DA ∆A
Ḃ = G(A, B) + DB ∆B.
A simple choice for the reaction part would be the activator-inhibitor model
from Gierer and Meinhardt, where species A is autocatalytic and activates B,
while B inhibits A. An even simplier choice is the activator-inhibitor model from
Schnackenberg, where the autocatalytic species A is the inhibitor of B and B
activates A. Both models form stripes in the Turing version and here we choose
the second one because it is mathematically easier to analyse:
F = k1 − k2 A + k3 A2 B
G = k4 − k3 A2 B.
k3 1/2 k3 1/2
u=A , v=B ,
k2 k2
DA t x DB
t= 2 , x= , d= ,
L L DA
k1 k3 1/2 k4 k3 1/2
a= , b= ,
k2 k2 k2 k2
L2 k2
γ=
DA
Note that the variables u and v are positive since they are concentrations of
reactants. By introducing the variables above, the system is described as follows
u̇ = γ(a − u + u2 v) + ∆u
| {z }
=:f (u,v)
v̇ = γ(b − u2 v) + d∆v
| {z }
=:f (u,v)
with the ratio of the diffusion d and the relative strength of the reaction versus
the diffusion terms γ which scales as γ ∼ L2 .
~ = (u−u∗ , v−
We start with a linear stability analysis of the reaction part using W
∗ ~ ∗ ∗
v ). We denote the steady state with W = (u , v ) and a partial derivative with
∗
fu = ∂f
∂u etc. This yields
~˙ = γAW
W ~
with the matrix !
fu fv
A= | .
gu gv W~ ∗
Linear stability is guaranteed if the real part of the eigenvalues λ is smaller than
zero, Re λ < 0. Thus, the trace of A is smaller than zero
tr A = fu + gv < 0 (8.122)
224
and the determinant larger than zero
We now turn to the full reaction-diffusion system and linearize it about the steady
state
~˙ = γAW
W ~ + D∆W
~
!
1 0
with D = .
0 d
In order to obtain an ODE from this PDE, we use the solutions of the Helmholtz
wave equation
∆W ~ + k2 W
~ =0
225
an instability, Re λ > 0, the constant part has to be negative. Since the first and
last terms are positive, this implies
This is the main result by Turing: An instability can occur if one component
diffuses faster than the other. (8.124) is only a necessary, but not a sufficient
condition. We require that the constant term as a function of k 2 has a negative
minimum.
(dfu + gv )2
> det A = fu gv − fv gu (8.125)
4d
The critical wavenumber can be calculated to be
1/2
det A
kc = γ
dc
with the critical diffusion constant from
for large t. Note however, that in this case also non-linear effects will become
important and thus will determine the final pattern.
In summary, we have found four conditions (8.122) - (8.125) for the Turing in-
stability to occur. We now analyze the Schnackenberg-model in one spatial di-
mension. We already noted that a + b > 0 and b > 0 for the steady state to make
sense. From the phase plane we see that f > 0 for large u and f < 0 for small u.
Hence, fu > 0 around the steady state. Thus, b > a.
From condition (8.122) and (8.124), we now calculate that d > 1 in this case (the
activation B diffuses faster in this model). In general, the conditions (8.122)-
(8.125) define a domain in (a, b, d)−space, the Turing space, in which the insta-
bilities occurs. The structure of the matrix A tells us how this will happen: as u
or v increases, u increases and v decreases. So, the two species will grow out of
phase. If there is a fluctuation to a larger A-concentration, it would grow due to
the autocatalytic feature. Locally this would inhibit B and it decreases strongly.
However, because B is diffusing fast, it now is depleted from the environment and
there A is not activated anymore. Therefore A goes down in the environment,
whereas B is high. This is the basic mechanism for stripe formation.
In the Gierer-Meinhardt model, the two species grow in phase. When the auto-
catalytic species A grows, so does B, because A is the activator in this model.
Now B diffuses out and suppressed A in the environment. This is an alternative
mechanism for stripe formation.
226
nπ
The domain size p has to be large enough for a wavenumber k = p to be within
the range of the unstable wavenumbers (γ ∼ L2 ):
2
nπ
γL(a, b, d) < < γM (a, b, d)
p
where L and M can be calculated exactly. Typically, the mode which grows has
n = 1. Whether the left or right solution occurs depends on the initial conditions.
If the domain grows to double size, than γ changes by four (γ ∼ L2 ). p stays
the same because it is measured in units of L. Now, the mode n = 2 behaves as
shown in the following figures. On this way, a growing system will develop a 1D
stripe pattern.
227