«I. INTRODUCTION AND SUMMARY Introductory courses in statistical physics often place their greatest emphasis on average quan- tities measured in ...»
Elementary simulation of tethered Brownian motion
John F. Beausang,1 Chiara Zurla,2 Luke Sullivan,3 Laura Finzi,2 and Philip C. Nelson1, ∗
Department of Physics and Astronomy, University of Pennsylvania, Philadelphia PA 19104, USA
Department of Physics, Emory University, Atlanta GA 30322
Ursinus College, Collegeville, PA 19426-1000
(Dated: October 1, 2006)
We describe a simple numerical simulation, suitable for an undergraduate project (or
graduate problem set), of the Brownian motion of a particle in a Hooke-law potential well.
Understanding this physical situation is a practical necessity in many experimental contexts, for instance in single molecule biophysics; and its simulation helps the student to appreciate the dynamical character of thermal equilibrium. We show that the simulation succeeds in capturing behavior seen in experimental data on tethered particle motion.
PACS numbers: 05.20.-y, 05.40.Jv, 87.15.Aa
I. INTRODUCTION AND SUMMARYIntroductory courses in statistical physics often place their greatest emphasis on average quan- tities measured in thermodynamic equilibrium. Indeed, the study of equilibrium gives us many powerful results without having to delve into too much technical detail. This simplicity stems in part from the fact that for such problems, we are not interested in time dependence (dynamics);
accordingly dissipation constants like friction and viscosity do not enter the formulas.
There are compelling reasons, however, to introduce students to genuinely dynamical aspects of thermal systems as early as possible—perhaps even before embarking on a detailed study of equilibrium1. One reason is that students can easily miss the crucial steps needed to go from a basic appreciation that “heat is motion” to understanding the Boltzmann distribution, and thus can end up with a blind spot in their understanding of the foundations of the subject. Although any kind of rigorous proof of this connection is out of place in a ﬁrst course, nevertheless a demonstration of how it works in a sample calculation can cement the connection.
A second reason to give extra attention to dynamical phenomena is the current surge in student interest in biological physics. Much current experimental work studies the molecular processes of life, or their analogs, at the single-molecule level, where simple mathematical descriptions do seem to capture the observed behavior.
One familiar setting where simple models describe statistical dynamics well is the theory of the random walk, and its relation to diﬀusion. Ref.1 showed an attempt to present classical statistical mechanics starting from the random walk, building on earlier cla
distribution must cancel changes caused by drift from a constant external force ﬁeld (gravity).
In this article we discuss a generalization of free Brownian motion that is important for interpreting a large class of current experiments in single-molecule biophysics: The Brownian motion of a particle tethered to a Hookean spring. Speciﬁcally, we investigate the dynamics of ﬂuctuations of such a particle in equilibrium. Although the theoretical analysis of this problem does appear in some undergraduate textbooks, it sometimes appears forbiddingly complex (e.g. Ref.3 ). We have found, however, that numerical simulation of the system is well within the range of an undergraduate project. Either as a project or a classroom demonstration, such a simulation brings insight into the emergence of equilibrium behavior from independent random steps, and also can serve as an entry into the topic of equilibrium ﬂuctuations.
The appendix gives an implementation written by one of us, an undergraduate at Ursinus College, in matlab.
II. EXPERIMENTAL BACKGROUND
As motivation, we brieﬂy mention two contexts in which Brownian motion in a harmonic (Hookelaw) trap has played a role in recent biological physics experiments.
Optical trapping1,4 is now an everyday tool for the manipulation of micrometer-scale objects (typically a polystyrene bead), and indirectly of nanometer-scale objects attached to them (typically DNA, RNA, or a protein). In this method, a tightly focused laser spot creates a restoring force tending to push a bead toward a particular point in space. When the trapping beam has a Gaussian proﬁle, the resulting force on the bead is often to a good approximation linear in bead displacement. Thus the bead executes Brownian motion in a harmonic potential well. In such a well the motions along the three principal axes of the well are independent.
The bead’s motion in one or two dimensions can be tracked to high precision, for example by using interferomotry, thus yielding a time series. The probability distribution of the observed bead locations then reﬂects a compromise between the restoring force, pushing the bead to the origin, and thermal motion, tending to randomize its location. The outcome of this compromise is a Gaussian distribution of positions, from which we can read the strength of the harmonic restoring force (“trap stiﬀness”). For practical reasons, however, it is often more accurate to obtain both trap stiﬀness and the bead’s eﬀective friction constant from the autocorrelation function of the bead position (see Sect. IV below). For example, slow microscope drift can spoil the observed probability distribution function (see Ref.4 ).
Our second example concerns tethered particle motion5. In this technique, a bead is physically attached to a “tether” consisting of a long single strand of DNA. The other end of the tether is anchored to a microscope slide and the resulting bead motion is observed. Changes in the bead’s motion then reﬂect conformational changes in the tether, for example the binding of proteins to the DNA or the formation of a long-lived looped state. Fig. 1a shows example data for a situation 3
FIG. 1: Left, sample experimental data for the x-component of the motion of a bead of radius Rbead ≈ 240 nm, attached to a DNA tether of length Ltether ≈ 3500 basepairs, or ≈ 1200 nm. The experiment, and the protocols used to remove drift from the raw data, are described in detail in Ref.6. For clarity only the ﬁrst 200 s of data are shown; the full data run lasted 600 s. Right, logarithm of the autocorrelation function of x expressed in nm2 (see text). Solid line, experimental data. Dots, simulation described in this paper, using parameters Aeﬀ = 72 nm and eﬀective viscosity 2.4 times that of water in bulk. [f:tpm]
FIG. 2: Left, histogram of measured bead x position for the experimental data shown in Fig. 1. The solid curve shows a Gaussian distribution with the same normalization and variance. Right, similar histogram for our simulation.
[f:histo] where such conformational changes are absent, that is, simple tethered particle motion.
As in the optical trap case, one can discard the dynamical information in the time series by making a histogram of particle locations. Fig. 2a shows the frequencies of occurence of various values of x. Rather detailed agreement between theory and experiment has been obtained for these histograms, including the slight deviation from Gaussian distribution shown in the ﬁgure6,7. In the present note, however, we are interested in a less sophisticated treatment of a more general question: Can we understand at least some aspects of the dynamical information contained in data like those in Fig. 1?
To this end, Fig. 1b shows the logarithm of the autocorrelation function, C(τ ) ≡ x(t + τ )x(t), where the brackets denote averaging over t. At τ = 0 this quantity is just the meansquare displacement, which would diverge for a free particle but instead has a ﬁnite value determined by the equipartition theorem. At large times the autocorrelation falls to zero, because two independent measurements of x are as likely to lie on opposite sides of the tethering point as they are to lie on the same side. In fact, the autocorrelation function should fall exponentially with τ, as we see it 4 does in Fig. 1b.3
We wish to simulate the motion of a bead of radius Rbead, attached to a tether of length Ltot, and compare our results to experimental data. For this we will need to know a speciﬁc property of DNA in typical solution: Its “persistence length” A is A ≈ 45 nm.1 We will suppose that the external forces acting on the bead are a hard-wall repulsion from the microscope slide, a tension force from the tether, and random collisions with surrounding water molecules (see Sect. V for further discussion). The tension force produces an eﬀective potential well that keeps the bead close to its attachment point. In fact, at low relative extension, the
tension exerted by a semiﬂexible polymer like DNA is approximately given by a Hooke-type law:1
f = −κx. The eﬀective spring constant is κ = 3kB T /(2ALtot ), where kB T ≈ 4.1 · 10−21 J is the thermal energy at room temperature and Ltot is the contour length of the polymer. Thus, again the motions in each of the x, y, and z directions are independent. Because the microscopy observes only x and y, we can reduce the problem to a two-dimensional one, and hence forget about the hard-wall force, which acts only along z. In fact, we will reduce it still further, by examining only the x coordinate of the bead position.
There is a subtlety in that we do not directly observe the endpoint of the polymer in an experiment. Rather, we observe the image of the bead; the image analysis software reports the location of the bead center, a distance Rbead from the attachment point. Thus the time series in Fig. 1a reﬂects the motion of the endpoint of a composite object, a semiﬂexible polymer attached by a ﬂexible link to a ﬁnal, stiﬀ segment of length Rbead. To deal simply with this complication, we note that a semiﬂexible polymer can also be approximately regarded, for the purposes of ﬁnding its force–extension relation, as a chain of stiﬀ segments of length 2A. In our case 2A ≈ 100 nm is not much larger than Rbead, so we will approximate the system as a single Hookean spring, with eﬀective length Ltot = Ltether + Rbead and an eﬀective persistence length Aeﬀ slightly larger than A. In the data we present, Ltether ≈ 3500 basepairs, or ≈ 1200 nm. We will ﬁt the data to obtain Aeﬀ.
Suppose we observe the bead at time intervals of ∆t. Without the tether, the bead would take independent random steps, each a displacement drawn from a Gaussian distribution with meansquare step length 2D(∆t), where D is the bead’s diﬀusion constant. If the bead were subjected to a constant force f (for example gravity), we could get its net motion by superimposing an additional deterministic drift on the random steps: ∆driftx = f /ζ. The friction constant ζ is related to D by the Einstein and Stokes relations: ζ = kB T /D = 6πηRbead, where η is the viscosity of water.
For the tethered case, at each step we instead use a position-dependent force −κx, where x is the current displacement. For small enough ∆t (perhaps smaller than the actual video frame rate), x will be roughly constant during the step, justifying this substitution.
5 Here again we ﬁnd a subtlety: The presence of the nearby wall creates additional hydrodynamic drag on the bead8,9. Moreover, the tether itself incurs signiﬁcant hydrodynamic drag impeding the system’s motion. Again for simplicity, we acknowledge these complications by ﬁtting to obtain an eﬀective viscosity ηeﬀ, which we expect to be larger than the value 10−3 Pa s appropriate to water in bulk.
IV. SIMULATION RESULTS[S:SR]
To summarize, the simulation implements a Markov process. Each step is the sum of a random, diﬀusive component with meansquare 2D(∆t), and a drift component −Dκx/kB T. The constants D = kB T /(6πηeﬀ Rbead ) and κ = 3kB T /(2Aeﬀ Ltot ) contain two unknown ﬁt parameters, the effective persistence length Aeﬀ and viscosity ηeﬀ. The output of the simulation is the probability distribution of positions, and the autocorrelation function, which may be compared to experimental data.
The simulation is deemed successful to the extent that the two ﬁt parameters take values reasonably close to the expected values, diﬀering in the expected directions, and the full functional forms of the outputs agree with experimental data. Figures Fig. 1b and Fig. 2 show that indeed the simple model works well. Our simulation took ∆t = 0.625 ms, for a total of about a million steps, which were then sampled every 40 ms to resemble the experimental data.
Our mathematical model made some naive simpliﬁcations. Two which have been mentioned involve the role of the bead radius, and certain sources of drag on the bead. In addition, there is a time scale for rearrangements of the DNA needed to change its extension, and for rotatory diﬀusion of the bead, which changes the location of the attachment point relative to its center. All of these eﬀects have been assumed to be lumped into eﬀective values of the ﬁt parameters.
Despite these simpliﬁcations, however, we did obtain two key qualitative aspects of the experimental data as outputs from the model, not set by hand. First, the autocorrelation function of equilibrium ﬂuctuations has the expected exponential form. Moreover, as a consistency check, the experimental histogram of bead positions had roughly the Gaussian form we would expect for motion in a Hookean potential well. Both of these results emerge as statistical properties of a large number of simple steps, each involving only a diﬀusive step combined with a drift step based on the current bead location.
The insights obtained from this exercise are diﬀerent from those obtained from the analytical solution; for example, students can see the average behavior emerging from the random noise as the simulation size grows. In addition, the simulation approach opens the door to replacing the assumption of a harmonic potential by any other functional form.
We thank I. Kulic and R. Phillips for many discussions.