[SOLVED] DNA matlab chain 4G1. MATHEMATICAL BIOLOGY OF THE CELL

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 = r(t) kr(t) + fr(t)
with E[fr(t)] = 0

E[fr,i(t1)fr,j(t2)] = 2kBTij(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 +

2Dtbi

where D = kBT

, ri = r(t = it) 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. Lets
call these , and respectively. You then define t = t/ the dimensionless time and ri = r

(
t =

it
)
/ the dimensionless position at the dimensionless discretized time it. 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:

ri+1 =
(

1t
)
ri +

2t 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 it 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 rn,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 N1:

rn,i+1 = rn,i t
[
k f(rn,i rn1,i) + k f(rn,i rn+1,i)

]
+

2t bn,i

r0,i+1 = r0,i t
[
k f(r0,i r1,i) + p

]
+

2t b0,i

rN,i+1 = rN,i t
[
k f(rN,i rN1,i) p

]
+

2t 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 = px and
calculate the chains extension e = E

[(
rN r0

)
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.
e, 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.