Brownian dynamics simulation of analytical ultracentrifugation experiments

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][3][4] and thematic issues of other journals [5][6][7].
In the AUC, particles move under influence of a centrifugal field, caused by rotation of the sample with angular velocity ω, which produces a centrifugal force (corrected by buoyancy) equal to  , 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), v is the partial specific volume of the solute particles and r 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: If this were the only action on the solute particles, their motion would be purely deterministic. If 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 Even if the initial (loading) concentration in the AUC cell is uniform, i.e., constant from the meniscus to the bottom in the AUC cell, which are placed at distances 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.
The AUC experiment yields the concentration profile as a function of time, i.e., the two-variables function 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.

The simulated system
We consider a solution of concentration 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 c = (r b − r m )/N r . Thus the cell is divided in slices of width c and transversal area jh, where j 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 jhcr 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 .
The mass concentration c(r i , t j ) in the real system is proportional to the number concentration in the simulated system: where Q is a proportionality constant which can be fixed imposing the condition for the initial uniform concentration over the whole volume of the cell, so that combining eq. 7 with 8 we readily find, 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
The main aspect of the Brownian dynamics (BD) simulation is the propagation of the particle trajectories. As we are considering non-interacting particles, the trajectory of one particle is independent of that any other, so that those quantities are the diffusion coefficient, 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 Although -as it will be verified later on -the end effects, at the solution meniscus and the bottom of the cell are of minor importance we take them into account. As the description sedimentation and Brownian motion near boundaries or walls seems problematic, we adopt 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 ) = 2r 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 During the large step the sedimentation velocity changes as r changes, but this change is deterministic, and as mentioned above the sedimentation drift is easily integrated as indicated in eq. 12 while, thanks to the fractal nature of the Brownian motion, the Brownian step follows the same law over the long time, Δ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. The initial position of the particle is assigned according to the uniform concentration in the sector-shaped cell. As the number of particles in a slice of thickness c 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 The initial position should be a random number obeying to the probability density expressed by eq. 14, and accumulative distribution function Thus we generate a random number with uniform distribution 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 j = 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.
Simulations were done for four different molecules covering an extremely large range of sedimentation to diffusion coefficients (i.e., in a very wide range of molecular mass), so that the extreme cases of vanishing and intense diffusion broadening are considered in our study. This includes a quite small molecule -g-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][14][15][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

Results and discussion
First of all, we checked that the results from the algorithm with a small number of large time steps, eq. 11 give practically the same results than that with many small steps, eq. 10, thus confirming the efficiency of the former. This was evident in all cases simulated, as shown for example in Figure 2.
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 Lammequation solver to determine the parameters by means of powerful optimization algorithms. In the singlespecies 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.
The simulated runs in the equilibrium mode were conducted for long times, following the evolution of the simulation profile to make sure that a further extension of the run time would not change the resulting final equilibrium 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 r m 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 Conditions for simulation of sedimentation equilibrium experiments, and molecular weight recovered from the analysis of the sedimentation profiles using eq. 17. Numbers between parenthesis are as in Table 1.

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 noninteracting 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/  Additional file 2: Sedimentation equilibrium experiment of cyclodextrin: Movie showing the trajectory followed by 1000 molecules of cyclodextrin, represented as red spheres in a typical ultracentrifuge cell, r m = 5.8 cm, r b = 7.2 cm in a sedimentation equilibrium experiment with ω = 30000 rpm and 300 000 s, at 20°C. The white spheres represent the water solvent.