# Brownian dynamics simulation of analytical ultracentrifugation experiments

- AI Díez
^{1}, - A Ortega
^{1}and - J Garcia de la Tore
^{1}Email author

**4**:6

**DOI: **10.1186/2046-1682-4-6

© Diez et al; licensee BioMed Central Ltd. 2011

**Received: **11 February 2011

**Accepted: **2 March 2011

**Published: **2 March 2011

## Abstract

### Background

We have devised a protocol for the Brownian dynamics simulation of an analytical ultracentrifugation experiment that allows for an accurate and efficient prediction of the time-dependent concentration profiles, *c*(*r, t*) in the ultracentrifuge cell. The procedure accounts for the back-diffusion, described as a Brownian motion that superimposes to the centrifugal drift, and considers the sector-shaped geometry of the cell and the boundaries imposed by the meniscus and bottom.

### Results

Simulations are carried out for four molecules covering a wide range of the ratio of sedimentation and diffusion coefficients. The evaluation is done by extracting the molecular parameters that were initially employed in the simulation by analyzing the profiles with an independent tool, the well-proved SEDFIT software. The code of simulation algorithm has been parallelized in order to take advantage of current multi-core computers.

### Conclusions

Our Brownian dynamics simulation procedure may be considered as an alternative to other predictors based in numerical solutions of the Lamm equation, and its efficiency could make it useful in the most relevant, inverse problem, which is that of extracting the molecular parameters from experimentally determined concentration profiles.

## Background

Since the invention of the analytical ultracentrifuge by Svedberg [1], the technique of analytical ultracentrifugation (AUC) has been a classical - and, thanks to advances in instrumentation and analysis software, it is still a most modern - technique for characterization of macromolecules and nanoparticles in solution. The reader may grasp the recent importance of this field in monographs [2–4] and thematic issues of other journals [5–7].

*ω*, which produces a centrifugal force (corrected by buoyancy) equal to ${\omega}^{2}rm(1-\overline{v}\rho )$, where

*r*is the instantaneous distance from the particle to the rotation axis,

*m*is the mass of the particle (

*m*=

*M*/

*N*

_{ A }where

*M*is the molecular weight and

*N*

_{ A }is Avogadro's number), $\overline{v}$ is the partial specific volume of the solute particles and

*ρ*is the solution density (nearly equal to the solvent density, if the solution is dilute). The velocity that the solute particles may acquire due to this effect is proportional to the centrifugal acceleration,

*υ*=

*sω*

^{2}

*r*, where the

*s*is the sedimentation coefficient, and modulated also by the friction coefficient

*f*of the particle in the viscous solvent:

*r*(

*t*) is the radial position of the particle at time

*t*, one easily finds (considering that

*υ*=

*dr*/

*dt*) that the position after some time, Δ

*t*, would be given by

*r*

_{ m }and

*r*

_{ b }, respectively, from the rotor, centrifugation will provoke some transport of the solute particles, and therefore a concentration gradient will be produced. This gradient will, in turn, generate a counterflow of solute in the direction of decreasing concentration, i.e., contrary to the centrifugal velocity. Macroscopically, at a point

*r*the counterflow would be determined by the first law of Fick,

*J*= −

*D*∇

*c*(

*r, t*), where

*D*is the diffusion coefficient, related to both

*f*and

*s*by the Einstein and Svedberg equations, respectively:

where *k*_{
B
} is Boltzmann constant, *T* is the absolute temperature and *R* = *k*_{
B
}*N*_{
A
} is the constant of perfect gases.

*c*(

*r, t*). It is from the analysis of this function that the information of interest about the solute particles, i.e., the values of

*s*,

*D*and

*M*can be determined. When the diffusional counterflow can be ignored - e.g., when the sedimentation is carried out at extremely large

*ω*, or for massive particles having a great

*s*/

*D*ratio - such analysis is simply based on eq. 2, as described in textbooks [8], and provides a value of

*s*. However, in general cases, particularly when the diffusional effect is influential, and therefore one could determine not only

*s*, but also

*D*, and

*M*therefrom, a rigorous consideration of both effects is required. The time-dependent concentration profile

*c*(

*r, t*) is determined by the balance between the centrifugal and diffusional effects. Classically, this has been expressed in macroscopic terms, in the form of the so-called Lamm equation [9], which expressed the balance between the centrifugal drift and the backwards flux that, according to the Fick law, has to occur because of the generated concentration gradient. The Lamm equation reads:

Eq. 6 is written in cylindrical coordinates, because the AUC geometry is radial, and sector-shaped cells are used (devised so that the radial trajectories would not collide with the lateral walls).

In spite of the basic nature of the concepts involved, the solution of the Lamm partial differential equation, with the geometry and boundary conditions of the AUC experiment is extremely difficult, and requires numerical methods [10, 11] based on time discretization and a finite-elements description of the AUC cell. [10]. Nonetheless, in modern AUC analysis programs, like SEDFIT [12], a Lamm-equation solver is embodied, enabling the prediction of computed *c*(*r, t*) for estimations of *s*, *D* and *M*, which are to be optimized as to fit the experimental *c*(*r, t*). In these procedures (apart from the fitting or optimization algorithms, a central piece is the prediction of *c*(*r, t*).

In this paper we investigate an alternative predictor of AUC concentration profiles. Instead of starting from the balance of the macroscopic flows established by the Lamm equation, we consider a microscopic description of the motion of particles under the simultaneous effect of a deterministic force, and the random forces characteristic of Brownian motion, so that the latter replaces the Fick's law description of macroscopic diffusion. In a simple (and somewhat naïve) approach, we formulate a Brownian dynamics algorithm to simulate the trajectories of particles. With our microscopic perspective, our method also discretizes time (Brownian simulation steps) and instead of finite elements use discrete particles. Carrying out such simulation for a sufficiently large number of particles, we can determine the time-dependent concentration profile in the AUC cell. In the next section we describe the procedure and demonstrate adequacy of its results, and finally we shall discuss on its performance, eventual advantages and possible extensions for further applications.

## Methods

### The simulated system

*c*

_{0}, initially uniform at

*t*= 0, in a sector-shaped cell (Figure 1). In the simulation, the solute is represented by a collection of a large number of particles

*N*

_{ part }. The radial distance is discretized in

*N*

_{ r }intervals of width

*χ*= (

*r*

_{ b }−

*r*

_{ m })/

*N*

_{ r }. Thus the cell is divided in slices of width

*χ*and transversal area

*ϕh*, where

*ϕ*is the angle of the sector comprising the cell, and

*h*its height, so that the volume of a slice placed at distance

*r*

_{ i },

*i*= 1, ...

*N*

_{ r }, is

*ϕhχr*

_{ i }. The concentration profile is to be calculated during a total time

*T*, at a series of

*N*

_{ t }time intervals, of duration

*τ*, so that

*T*=

*N*

_{ t }

*τ*and

*t*

_{ j }=

*jτ*,

*j*= 1, ...

*N*

_{ t }.

*c*(

*r*

_{ i },

*t*

_{ j }) in the real system is proportional to the number concentration in the simulated system:

*Q*is a proportionality constant which can be fixed imposing the condition for the initial uniform concentration over the whole volume of the cell, $\frac{1}{2}\varphi h({r}_{b}^{2}-{r}_{m}^{2})$:

Note that, apart from several constants relative to the instrument or the simulation the concentration is not determined by *n*(*i, j*), but by *n*(*i, j*)/*r*_{
i
}. This takes into account what is called in the AUC terminology, the radial dilution effect.

What we simulate corresponds to a highly diluted solution, in which there are no particle interactions. Thus we can generate trajectories of individual particles, independent of each other.

### Brownian dynamics simulation

*D*, and the force

*F*, acting on the particle. In the present problem, we could begin with a basic algorithmic procedure in which the time step,

*δt*is rather small, so that neither of the quantities determining the step change appreciably during

*δt*. Then, the running algorithm for the particle's position would be:

where *δr*_{
sed
} = *ω*^{2}*srδt* is the deterministic sedimentation drift of the particle with instantaneous, position-dependent velocity *ω*^{2}*sr*, while the random Brownian displacement has zero mean and variance $<\delta {r}_{Brow}^{2}>=2D\delta t$.

*ad hoc*criteria. As for the meniscus, if after the step

*r*<

*r*

_{ m }, we set

*r*=

*r*

_{ m }. Regarding the bottom, if

*r*>

*r*

_{ b }, the particle had hit the bottom of the cell during the step; then we assume it should bounce and correct the position, taking

*r*− 2(

*r*−

*r*

_{ b }) = 2

*r*

_{ b }−

*r*. After testing that this algorithm, in which the trajectory is divided in a very large number of small time steps, predicts correctly the concentration profiles (see below) we intended to devise a procedure with larger times steps, which would be computationally faster. The displacement over a large time step Δ

*t*is the result of the integration of the small increments in eq. 10, so we can write

*r*changes, but this change is deterministic, and as mentioned above the sedimentation drift is easily integrated as indicated in eq. 12

*r*

_{ Brow }being a random number of zero mean and variance

Thus the algorithm based on eqs. 11, 12 and 13 could be applicable to arbitrarily large time steps (even as large as the time interval *τ* between registers). This is essentially true if there were no end effects, i.e., in infinite, unbound AUC cell. For the sake of simplicity, we still adopt the simple criteria that particles stop at the meniscus and bounce at the bottom. Thus the only defect introduced by this procedure would be an inaccurate prediction of the concentration near the meniscus and bottom. In this regard, we note that the end-effects also affect other prediction procedures, like those based in Lamm-equation solvers, and influence the experiment itself, so that it is a common practice to discard the two terminal regions in the analysis of AUC experiments.

### Procedure

Summarizing from the previous description, Brownian dynamics trajectories are simulated for a large number of particles, *N*_{
part
}. The trajectory of one particle is monitored, determining at successive times *t*_{
j
} the interval of radial position *r*_{
j
}. Then the counter for those interval and position is increased *n*(*i, j*) → n(*i, j*) + 1.

*χ*is proportional to r, the probability of having (in the uniform solution) a particle at a distance

*r*is

*p*(

*r*) ∝

*r*. Integrating from

*r*

_{ m }to

*r*

_{ b }gives

*u*∈ (0, 1) which is equated to

*F*(

*r*) to obtain the initial position.

If *t*_{
run
} is the total time of the experiment (usually, several hours), the simulation for each particle consists of a number of steps *N*_{
steps
} = *t*_{
run
}/Δ*t*. For recording purposes, a number of registers (scans in the AUC experiment), equal to *N*_{
scans
}*T* = *t*_{
run
}/*τ* is made, observing the radial distance, and therefrom the index *i*, of the slice at which the particle is at that instant, *j* = *t*_{
run
}/*τ*. Then one unit is added to the counter, *n*(*i*, *j*). Trajectories are simulated for a sufficiently large number of particles *N*_{
part
}. At the end of the simulations, concentration profiles *c*(*r*_{
i
}*, t*_{
j
}) are obtained from the *n*(*i, j*) counter using eq. 9.

### Simulation data

The data employed in the simulation reflect those of real AUC experiments. The geometry of the cell is given by *r*_{
m
} = 5.80 cm, *r*_{
b
} = 7.20 cm and sector angle *ϕ* = 3 degrees (0.05 radians). The rotor speed *ω* and the duration of the experiment *t*_{
run
} was varied depending on the molecule being simulated, and the mode (velocity or equilibrium) of the experiment. In some cases also the meniscus position was varied for equilibrium experiments. We found that *N*_{
part
} = 10^{5} particles suffices to obtain a rather low level of noise as it will be shown below.

*γ*-cyclodextrin - a small globular protein - lysozyme - a moderately long flexible polymer - poly(ethylene oxide) - and a very long and stiff DNA from T7 bacteriophage. Values for

*s*,

*D*and

*M*are taken from the literature [8, 13–16], and listed in Tables 1 and 2, along with the conditions regarding spinning time and velocity for each sample in both modes, for velocity and equilibrium sedimentation experiments

Experimental properties

Molecule | γ-cyclodextrin | Lysozyme | Polyoxyethylene | DNA T7 |
---|---|---|---|---|

$1-\overline{v}\rho $ | 0.333 | 0.298 | 0.170 | 0.450 |

10 | 1.50 | 14.3 | 300 | 26 000 |

| 0.53 | 1.89 | 2.30 | 31.8 |

10 | 27.0 | 10.8 | 1.10 | 0.066 |

| 55 000 | 40 000 | 40 000 | 10 000 |

| 40 000 | 40 000 | 40 000 | 60 000 |

10 | 1.44 (-4) | 13.9 (-3) | 300 (0) | 23 278 (-10) |

| 0.54 (+2) | 1.84 (-3) | 2.26 (-2) | 31.4 (-1) |

10 | 27.4 (+1) | 10.8 (0) | 1.08 (-2) | 0.073 (+11) |

Equilibrium experiments

Molecule | γ-cyclodextrin | Lysozyme | Polyoxyethylene |
---|---|---|---|

| 30 000 | 15 000 | 5 000 |

| 300 000 | 500 000 | 500 000 |

10 | 1.27 (-16) | 12.6(-12) | 250.0 (-17) |

## Results and discussion

In order to evaluate the accuracy of the results for the prediction of sedimentation profiles in the velocity mode, we have compared them with calculations with the well-known SEDFIT software [12]. This tool, which includes a sophisticated Lamm-equation solver has an utility for predicting the concentration profiles. Comparison of our simulated profiles with the SEDFIT predictions resulted in a good agreement, as shown in Figure 2. The determination of concentration profiles has the final purpose of determining the molecular parameters of interest, so the most relevant evaluation of their prediction is the confirmation of whether their analysis provides correct values for those parameters. In order to do so with our simulation results, we employed the analysis tool of SEDFIT as an independent and robust criterion. Indeed, SEDFIT uses a reliable Lamm-equation solver to determine the parameters by means of powerful optimization algorithms. In the single-species mode of SEDFIT, the program provides the values of the molecular weight, *M*, and the sedimentation and diffusion coefficients, *s* and *D* of the sample. We made the SEDFIT analysis with typically 50 scans (*t* values) each covering 1400 radial positions *r* values.

In Table 1 we report the *s*, *M*, *D* and values for the four samples considered in our study. We note the high accuracy of the recovered values of these three parameters, reflected in their very small deviations (particularly in the case of *s*) from the values employed as input in the simulation. The agreement is good for the four cases, which, as indicated above, cover a wide range of the sedimentation-to-diffusion ratio and molecular mass, including in the velocity experiments a molecule as small as cyclodextrin.

*c*(

*r*) profile. An example is presented in Figure 3. For the analysis of the simulated runs in the equilibrium mode, we adopted the classical equation [8]:

which, as usual, is linearized in the form of a plot of ln *c*(*r*) vs. ${r}^{2}-{r}_{m}^{2}$ whose slope is ${\omega}^{2}M(1-\overline{v}\rho )/(2RT)$, from which the molecular weight *M* is extracted. The results are listed in Table 2. Here we observe that the recovered values of *M* are not as good as those obtained in the velocity mode, with deviations of about 15%. Among other reasons (like the smaller amount of information resulting from equilibrium experiments, and the errors inherent to the determination of the slope in the fit to the linearized equation) this may be because the simple eq. 17 neglects radial dilution.

## Computing details

The extreme simplicity of the algorithm that simulates the trajectory of one molecule makes the simulation scheme very well adapted for parallelization, thus taking benefit of present multi-core platforms, because each trajectory can be generated in a separate thread/core. We have implemented OpenMP directives in our Fortran 90 code and tested the performance in a DELL T5500 workstation with two Intel Xeon X5660 processors. 1000 simultaneous simulation with *N*_{
part
} = 10^{5}, *N*_{
r
} = 100 radial positions and *N*_{
t
} = 50 recorded scans took 46.5 CPU seconds. Broadly speaking, our algorithm, which is easily parallelized, is able to run one thousand *c*(*r, t*) calculations in less than one minute. The computing speed of the algorithm is crucial in its main use, namely, in the analysis of *c*(*r, t*) experimental profiles by any kind of fitting to computed profiles, and the speed of our procedure seems suitable for that purpose.

A nice feature of the Brownian simulation of ultracentrifugation is that it allows to visualize the simulated trajectories using computer graphics. This may be of utility for demonstrative purposes (e.g. in teaching AUC principles). We have produced two videos, showing the evolution of the solute particles in the AUC as the sedimentation proceeds, one at high rotor speed, in the velocity mode, and another at low speed, that reaches the sedimentation-diffusion equilibrium. These videos accompany this paper as additional files 1 and 2.

## Conclusions

In this work we have presented a simple procedure, based on a Brownian dynamics, microscopic simulation, for predicting the time/position-dependence of concentration (concentration profiles *c*(*r, t*)) during the AUC experiment, which can be regarded as an alternative to the numerical solution of the macroscopic Lamm equation. The correctness of the procedure has been tested comparing the profiles with those computed by the Lamm-equation solution, and the concordance of the molecular parameters recovered from them with those used in the simulation.

Having presented this proof-of-concept, some advantages of our scheme can be hinted - although they all remain to be evaluated in future work. An important aspect to be considered is computational efficiency. As commented above, our computational procedure has the feature of being perfectly parallelizable. A comparison with the numerical solution of Lamm solution requires further the labor, implementing codes for those procedures and making a side-to-side analysis of computing speed and requirements. Further work is planned in this regard.

Apart from computing efficiency, the BD scheme has the potential advantage that it can treat easily cases of arbitrary complexity. In our proof-of-concept we have restricted ourselves to the simplest case of identical, non interacting particles (extension to different but still non-interacting particles is trivial). An extraordinary utility of AUC is the characterization of macromolecular interactions. For problems with interacting particles, the BD scheme can be easily adapted. In BD, detailed interactions between the particles can be modeled whereas the continuum method only allows for averaged density dependent potentials. Even if computing efficiency would suffer in such, more complex problems, still the BD approach, which is based on first principles and allows explicit description of interaction between particles (or the effect of special conditions in the AUC experiments) may be a valuable tool for testing other approaches.

## Computer programs

Fortran 90 source code with OpenMP directives will be freely downloaded from our web site, http://leonardo.inf.um.es/macromol/

## Declarations

### Acknowledgements

This work was performed within a *Grupo de Excelencia de la Región de Murcia* (grant 04531/GERM/06). Support also provided by grant CTQ-2009-08030 from *Ministerio de Educación y Ciencia* (MEC), including FEDER funds. A.O. acknowledges a postdoctoral fellowship from *Universidad de Murcia*, and A.I.D. is recipient of a predoctoral fellowship from MICINN.

## Authors’ Affiliations

## References

- Svedberg T, Pedersen K: The Ultracentrifuge. 1940, Oxford: Clarendon PressGoogle Scholar
- Mächtle W, Börger L: Analytical Ultracentrifugation of Polymers and nanoparticles. 2006, Berlin, Heidelberg: Springer-VerlagGoogle Scholar
- Harding S, Rowe A, Horton J: Analytical Ultracentrifugation: Techniques and Methods. 1992, Cambridge, Chester: Royal Society of ChemistryGoogle Scholar
- Scott D, Harding S, Rowe A: Analytical Ultracentrifugation: Techniques and Methods. 2005, Cambridge: Royal Society of ChemistryView ArticleGoogle Scholar
- Solovyova A: 17th international symposium on analytical ultracentrifugation and hydrodynamics. Eur Biophys J. 2010, 39: 345-346. 10.1007/s00249-009-0406-4.View ArticleGoogle Scholar
- Cölfen H: Analytical Ultracentrifugation. Macromol Biosci. 2010, 10: 687-688.View ArticleGoogle Scholar
- Zhao J, Schuck P: Methods.Google Scholar
- van Holde KE, Johnson WC, Ho PS: Physical Biochemistry. 1998, New Jersey: Prentice HallGoogle Scholar
- Lamm O: Die Differentialgleichung der Ultrazentrifugierung. Ark Mat Astr Fys. 1929, 21B: 1-4.Google Scholar
- Claverie JM, Dreux H, Cohen R: Sedimentation of generalized systems of interacting particles. I. Solution of systems of complete Lamm equations. Biopolymers. 1975, 14: 1685-1700. 10.1002/bip.1975.360140811.View ArticleGoogle Scholar
- Dishon M, Weiss G, Yphantis D: Numerical solutions of the Lamm equation. I. Numerical procedure. Biopolymers. 1966, 4: 449-455. 10.1002/bip.1966.360040406.View ArticleGoogle Scholar
- Schuck P: Size distribution analysis of macromolecules by sedimentation velocity ultracentrifugation and Lamm equation modeling. Biophys J. 2000, 78: 1606-1619. 10.1016/S0006-3495(00)76713-0.View ArticleGoogle Scholar
- Pavlov G, Korneeva E, Smolina N, Shubert U: Hydrodynamic properties of cyclodextrin molecules in dilute solutions. Eur Biophys J. 2010, 39: 371-379. 10.1007/s00249-008-0394-9.View ArticleGoogle Scholar
- Müller J: Prediction of the rotational diffusion behavior of biopolymers on the basis of their solution or crystal structure. Biopolymers. 1991, 31: 149-160.View ArticleGoogle Scholar
- Luo Z, Zhang G: Scaling for Sedimentation and Diffusion of Poly(ethylene glycol) in Water. J Phys Chem. 2009, 113: 12462-12465.View ArticleGoogle Scholar
- Crothers D, Zimm B: Viscosity and Sedimentation of DNA from Bacteriophages T2 and T7 and Relation to Molecular Weight. J Mol Biol. 1965, 12: 525-536. 10.1016/S0022-2836(65)80310-2.View ArticleGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.