4G1. MATHEMATICAL BIOLOGY OF THE CELL
COURSEWORK 2
Write a report that addresses all the questions below. It should include figures and snapshots
of your simulations. Provide your codes in an appendix. The report should not exceed 5 pages
(excl. appendix) and is due Thursday, March 23 by 4pm on moodle. Submit a single PDF,
combining the coversheet and the report.
I. Brownian particle in a harmonic potential. The equation governing the dynamics of an
inertialess Brownian particle in a harmonic potential at time t is
(1) 0 = −ζ ṙ(t)− kr(t) + fr(t)
with E[fr(t)] = 0
E[fr,i(t1)fr,j(t2)] = 2ζkBTδijδ(t1 − t2)
where δij = 1 if i = j (0 otherwise) and i, j = x, y, z.
I.1. (5 points) Solve for E[r(t)2] with the initial condition r(0) = 0. Define τ = ζ/k.
I.2. (5 points) Discuss the limits of E[r(t)2] for t� τ and t� τ .
I.3. (15 points) Perform a Brownian Dynamics simulation to check your results.
Hint 1: The Matlab code BD_free_motion.m below computes the trajectory of a free particle
(k = 0) using an explicit Euler method with time step ∆t:
ri+1 = ri +
√
2D∆tbi
where D = kBT
ζ
, ri = r(t = i∆t) and bi is a random vector where each coordinate is independent,
normally distributed with zero mean and unit variance.
1 function result=BD free motion(dim,N steps) % dim: dimension
2 % N steps: number of time steps
3
4 delta t = 0.001; % Time step
5 diff = 1.0000; % Diffusion coefficient
6 Std BF = sqrt(2*diff*delta t); % Standard deviation of Brownian force
7
8 r = zeros(dim,N steps+1); % Initial array of positions
9
10 % Brownian Dynamics
11 for i=1:N steps
12 r(:,i+1) = r(:,i) + Std BF * normrnd(0,1,[dim,1]);
13 end
14
15 result = r; % the result
16 end
Hint 2: In fact, a computer can only deal with dimensionless numbers. For the Brownian
Dynamics simulation, you need to choose characteristic scales of time, length and energy. Let’s
call these θ, λ and � respectively. You then define t̃ = t/θ the dimensionless time and r̃i = r
(
t =
iθ∆t̃
)
/λ the dimensionless position at the dimensionless discretized time i∆t̃. A natural choice
1
2 COURSEWORK 2
for these scales in the case of the trapped particle described by Eq. (1) is θ = τ , λ =
√
kBT/k
and � = kBT . In that case, you may verify that the explicit Euler method reduces to:
r̃i+1 =
(
1−∆t̃
)
r̃i +
√
2∆t̃ bi
and find the theoretical dimensionless expression E
[(
r
(
t = θt̃
)
/λ
)2]
= E[r̃
(
t̃
)2
] using your result
from question I.1.
II. Force-extension of DNA. We want to study the mechanical response of a single DNA
molecule. We model the DNA as a bead-spring chain, made of N + 1 particles connected by
N springs. The springs all have a stiffness κ and length at rest `. The position of the bead
n ∈ {0 . . . N} at the discretized time i∆t is written rn,i.
p p
II.1. (10 points) It is now more convenient to use λ = `, θ = ζ`
2
kBT
, � = kBT as the characteristic
scales for the simulation (see Hint 2 above). We thus introduce r̃n,i = rn,i/λ and t = θt̃,
as well as the dimensionless spring constant k = κ`
2
kBT
. For the bead-spring model, show
that the dimensionless Euler integration method is written, with i ≥ 0 and 1 ≤ n ≤ N−1:
r̃n,i+1 = r̃n,i − ∆t̃
[
k f(r̃n,i − r̃n−1,i) + k f(r̃n,i − r̃n+1,i)
]
+
√
2∆t̃ bn,i
r̃0,i+1 = r̃0,i − ∆t̃
[
k f(r̃0,i − r̃1,i) + p̃
]
+
√
2∆t̃ b0,i
r̃N,i+1 = r̃N,i− ∆t̃
[
k f(r̃N,i − r̃N−1,i)− p̃
]
+
√
2∆t̃ bN,i
where f(a) = a(1− |a|−1) and p̃ = p`
kBT
is the dimensionless pulling force.
II.2. (40 points) Perform the simulation for various values of the pulling force p̃ = p̃x̂ and
calculate the chain’s extension ẽ = E
[(
r̃N − r̃0
)
· x̂
]
in the direction of the force (x̂ being
the unit vector in this direction). Show in particular how the force-extension curve p̃ vs.
ẽ, depends on k and N .
II.3. (15 points) Prove that the expression for the force-extension curve of the freely-jointed
chain with links of length ` is:
e
N`
= coth
(
p`
kBT
)
−
kBT
p`
Compare it with your simulation results and discuss the similarities and differences.
II.4. (10 points) You will find on moodle some experimental data for the force-extension of
λ-phage dsDNA, which has a contour length of L = 32.7 µm and a persistence length of
53 nm (i.e. a Kuhn length ` = 106 nm). Compare the data with the freely-jointed chain
model, and with the following interpolation for the worm-like chain model:
p`
2kBT
=
1
4
(
1−
e
L
)−2
−
1
4
+
e
L
In both cases, discuss the similarities and differences.