Note, as of 2007,a significant update and improvement in the Lamm equation solution is implemented, as described in

**P.H. Brown, P. Schuck (2007) A new adaptive
grid-size algorithm for the simulation of sedimentation velocity profiles in
analytical ultracentrifugation. **

**Click on this link to read how this works: **
Computer Physics Communications in
press

**This work is not yet included in this introductory page.**

The finite difference approach and its limitation

The finite element solution by Claverie

The moving frame of reference finite element solution

The Lamm equation describes the evolution of the concentration distribution
c(r,t) of a species with sedimentation coefficient s and diffusion coefficient D
in a sector-shaped volume and in the centrifugal field w^{2}r.
It is a partial differential equation that, in the simplest case, takes the
following form,

(1)

which can be derived from a mass conservation equation with the local
diffusion fluxes D(dc/dr) and the local sedimentation fluxes csw^{2}r.

Typical shapes of families of concentration distributions described by the Lamm equation look like this (centrifugal force pointing to the right, details would depend on the time-scale and rotor speed):

(a)(b)

(c)(d)

(e)(f)

[for species with initially uniform distribution, sedimenting at different values of s and D: (a) intermediate s and D; (b) intermediate s, low D; (c) high s, low D; (d) low s, high D; (e) flotation with intermediate s and D; (f) sedimentation with intermediate s and D with initial lamellar distribution of c(r)]

Some general features can be observed:

(1) the centrifugal field causes concentration gradients (d) or sedimentation boundaries (a, b, c, e);

(2) diffusion broadens this boundary with time (visible best in a, b, e, f);

(3) radial dilution (visible best in a, b, c) or radial concentration (for floatation in e), which is caused by the sector-shaped solution column;

(4) end effects of the solution column due to accumulation of the material at one end (omitted in a, b, c, but visible in d, e) and depletion from the other end of the solution column;

(5) after a long time, an exponential equilibrium distribution is approached (visible best at low speeds or high D, such as in d)

(6) it should be noted that, although the Lamm resembles a diffusion
equation, it has an extra transport term sw^{2}r
which is position dependent. It varies by about 15 - 20% for typical
experimental situations.

These features make the Lamm equation difficult to solve, and eliminate the possibility of a closed analytical solution. Since it was introduced in 1929, because of the great practical importance, several approximate analytical solutions have been derived for different limiting cases (see the monograph by Fujita, 1962).

The simplest approximation is that by Faxén (1929), which can be written as

(2)

with the meniscus position *r _{m}*, the boundary position of a
non-diffusing species

Improved analytical solutions have been developed, including higher-order terms of analytical serial expansions, which give excellent approximations for ranges of parameter values (s, D, w) resembling many commonly encountered experimental configurations (Holladay, 1979; Behlke, 1997; Philo, 1997).

A completely different strategy for solving the Lamm equation are numerical simulations. Instead of formulating approximate analytical expressions for the concentration distribution as a function of time, the numerical simulations aim to develop computer algorithms that mimic the evolution of the concentration distribution with time. Given any distribution c(r,t), the elementary step in the numerical algorithms is the calculation of the distribution at a later time c(r,t+Dt). Many different strategies have been developed (Bethune and Kegeles, 1961, Cann and Kegeles, 1974, Claverie, et al., 1975, Cox, 1965, Dishon, et al., 1966, Marque, 1992, Minton, 1992, Sartory, et al., 1976, Schuck et al, 1998, and Schuck, 1998).

In Sedfit, numerical simulations are implemented for the following reasons:

* Although the numerical simulations may require hundreds or thousands of iteration steps, these steps can be made computationally very efficient, so that there is no relevant difference in computation speed to the approximate analytical solutions.

* The approximate analytical expressions are based on truncated serial expansions that converge best for certain ranges of parameters (s, D, w, t) or initial conditions. In fact, many of the analytical solutions are based on idealizations such as infinite solution columns. In contrast, the numerical solutions, if properly designed, are naturally very flexible to adapt much better to different parameter values and initial conditions. For example, a well-designed algorithm can give equal precision for all curves (a) to (e) shown above, as well as many other situations. It is even possible to take experimental scans as initial conditions for the simulation of the evolution in time.

* The precision is easily controlled to be well below experimental errors of data acquisition for all cases. This can be achieved by sufficient discretization in the spatial variable and adaptive step-size control in the time variable.

* Numerical solutions allow the acceleration phase of the rotor to be included in the simulation.

* They can be naturally extended to take chemical reactions and concentration-dependent sedimentation into account.

However, numerical solutions can also have drawbacks: If they are not properly designed, they may converge only very slowly to the correct solution, and exhibit spatial and temporal oscillations in the concentration distributions. Therefore, in the following the different strategies and their advantages and disadvantages are briefly described.

**The finite difference approach and its limitation**

The simplest way of simulating Lamm equation solutions is the division of the radial variable into a series of compartments (i = 1...N, for example with equal width Dr). Consideration of the diffusion and sedimentation fluxes between the compartments allows to predict the evolution of the concentration in each compartment.

Example for the division of a boundary (red) into a series of compartments with different concentrations (visualized as rectangular boxes of different heights).

This is inspired by compartment models as commonly used in models for chemical kinetics, and the intuitive application of compartments for the spatial discretization.

If one writes the concentrations in each compartment in vector form c(r,t), the concentration changes occurring within a small time-interval Dt can be formally expressed as

(3a) , or (3b)

(with the identity matrix **I **and a propagation matrix **E**; this
follows from the linearity of the Lamm equation.) If we count the fluxes going into
and out of any compartment, we obtain for the discretized net diffusion fluxes -(D/Dr)[c(i+1)
- 2 c(i) + c(i-1)] and for the net sedimentation fluxes w^{2}[c(i-1)r(i-1)
- c(i)r(i)]. As a consequence, one arrives at the following form of the propagation matrix E

(4)

This is a tridiagonal system, and the evolution Eq. 3 can easily be solved for series of time-steps (Schuck et al., 1998).

One problem commonly encountered with such propagation schemes as in Eq. 3 is an oscillation caused by sequential too high and too low concentration estimates in neighboring compartments, which alternate in sign with each time-step. Once formed, such oscillations can often amplify. They can originate, for example, in regions of high concentration gradients.

A very effective way around this problem was described by Crank and Nicholson
(1947). The idea is to dampen the oscillations and get higher
precision of the time-steps by evaluating concentrations in the right-hand side
of Eq. 3a in the *middle during the
time-step D t*, approximated as (c(t)
+ c(t+Dt))/2.
The physical basis for this is that as the concentrations change gradually
during the time-step Dt, so do the fluxes, which are
proportional to the concentrations. Although we do not calculate the exact
form of these changes during the time-step Dt, it is
a good assumption that the average flux during Dt is
proportional to the average concentration during Dt,
which we approximate by the arithmetic mean of c(t) and c(t+Dt).
At first, it may seem that this is a circular argument, since we don't know c(t+Dt)
yet, but a simple inspection of the equation 3a shows that in fact (c(t)
+ c(t+Dt))/2 can be
used in place of c(t) in the right-hand side,
and the resulting equation can be easily rearranged to solve for c(t+Dt)
as

(5)

The Crank-Nicholson scheme gives a significant improvement in speed and accuracy of the simulation.

The finite difference approach is fast and works very well for small molecules (high D, or low
rotor speeds). However, for larger s-values with smaller diffusion
coefficients the approximation becomes poor, because the discretization into
compartments creates an artificial diffusion term of a magnitude proportional to
w^{2}sDr , where Dr
is the size of a single compartment (Schuck et al., 1998). For high D, the
extra diffusion term can be neglected, but for small D it cannot, unless a very
fine grid (small Dr) is used, which is very
inefficient.

The origin of this error can be seen as follows: Assume the ideal linear translation of a boundary in the absence of diffusion and radial dilution, as approximated by the compartment model.

At time 0 it may look like a perfect rectangular boundary, with empty compartments #1 and #2, and equally filled compartments #3 to #6.

After a time-increment sufficient for the boundary to move just half the size of the compartment, Dr/2, the concentration distribution would move as indicated by the blue dashed lines. However, in our discretization scheme this means that half the concentration of each compartment flows into the right-hand neighbor compartment, and based on this we calculate the new concentrations in each compartment. Because just half the concentration is left in compartment #3, it is half the height.

After another time-increment, again half of the concentration in each compartment flows to the right - leaving 25% in #3, 75% in #4, and 100% in #5 and #6.

It can be seen that instead of translating the rectangular boundary, the
compartment discretization acts like a *diffusion *of the boundary; the
interpretation of compartments as spatial elements introduces bias in the
spatial distribution.

This effect is well-known in numerical mathematics and referred to as 'numerical diffusion'. In the appropriate text-books, you can find that in order to avoid this 'numerical diffusion' the simulation may require several orders of magnitude finer grid size than the physical scales of primary interest (Durran, 1999).

In summary, the finite difference scheme is a good method for the Lamm equation only for small s and large D. The inherent problem is that the sedimentation term is not well described by the compartment model, which leads to poor convergence of this algorithm to the correct distribution, and requires very fine discretization. Although successive grid refinements can be incorporated into nonlinear regression, the finite element methods described in the following are significantly more efficient (Schuck et al., 1998).

**The finite element solution by Claverie**

Claverie has introduced the finite element technique to the problem of solving the Lamm equation, which is less empirically (or physically) motivated, and instead a more mathematically based approach. Nevertheless, it has many intuitive features. Instead of the approximation of the boundary by 'rectangular' compartments, the boundary is described as the mathematical sum of triangular-shaped elements called hat-functions:

(6)

where each of the hat-function looks like, for example, the red one indicated below (triangle in left figure).

As can be easily verified, they have the property that when neighboring hat-functions added up, they describe a step-wise linear function (e.g. the red line). Eq. 6 separates the spatial and the time variable, and the problem reduces to the question of how much of each element is required to describe the boundary well at a later time (such as the reduced magnitude of the red triangle in the right figure above).

The trick of getting equations for the amplitudes c_{k}(t) of each
triangle is the
multiplication of the Lamm equation with one element P_{k}(r) followed
by integration over rdr from meniscus to bottom (equivalent to the volume
integral in the sector-shaped cell). Integration by parts of the rhs can
be performed, exploiting the vanishing total flux at the meniscus and bottom (Claverie,
et al., 1975), and we arrive at

(7) [for all k]

Finally, we insert our approximation Eq. 6 (renaming our running variable to j) into Eq. 7, which leads to an equation system

(8a) [for all k]

It should be noted that since the hat-functions P_{k}(r) are simple,
the derivatives dP_{k}(r)/dr are very simple, too. Also, since the
spatial variable only occurs in the form of the hat functions, we
can easily calculate analytically the integrals for all (j,k), and abbreviate Eq.
8a to

(8b)

This is a linear equation system that can be easily solved. The matrices **A**^{(1)}, **A**^{(2)} and **B** are
tridiagonal (because the pair-wise integral over the hat-functions is only
non-zero for identical or neighboring hats), and their numerical values can be found here
(as C-code). Discretization of the time-variable and writing the concentrations
in vector-form leads to the equation system

(9a)

or in Crank-Nicholson scheme

(9b)

This is a triangular equation system very similar to the one derived above for finite difference system, but because of the use of hat-functions instead of the compartments, it does not introduce the artificial diffusion term!

For calculating the evolution in time, Eq. 9b can be used, as long as the time-steps Dt are not too large. Initially, or when rapid changes in the concentration profiles take place, a small Dt has to be used, while near equilibrium very large Dt may be tolerable (and preferable). Therefore, to keep the discretization error small but the simulation efficient, it is very useful to adjust the size of the time-step dynamically after each step (Schuck et al., 1998).

This simulation of sedimentation works extremely well for small and
medium-sized molecules; it provides a very high precision at only moderate
computational cost. The algorithm has problems, however, in the case of
large molecules, where the sedimentation flux is much larger than
diffusion. In this case, the sedimentation boundary is very steep, and
very small time-steps are required in order avoid too large changes in the c_{k}(t)
values when the boundary moves across the compartments. In fact, the
algorithm is unstable for the limiting case of no diffusion (where
even very small time-steps can lead to negative concentrations c_{k}(t),
followed by large oscillations). This is not very appealing, especially since there is a simple analytical solution of the Lamm equation for particles without
diffusion:

(10)

describing a rectangular step at the position *r _{m}exp(*w

**The moving frame of reference finite element solution**

The basic idea for an algorithm that solves the Lamm equation in a stable and efficient manner for all limiting cases is the movement of all radial grid (reference) points exactly like a non-diffusing particle (Schuck, 1998). Viewed from this reference frame, there is no sedimentation visible, just diffusion, which by itself is a comparatively simple computational problem. This is similar to the idea of transformations of the radial variable that was employed for the derivation of the approximate analytical solutions of the Lamm equation (Faxén, 1929 and others, see Fujita 1962), but incorporated here in the finite element framework.

As indicated in this schematic view, in this algorithm the hat functions at
the different grid points *sediment* with the particle, such that the
relative changes in the amplitudes of our hat functions remain small, which allows higher precision and efficiency at
large sedimentation flux. Mathematically, we can formulate this idea by
approximating the
sedimentation profiles with a linear combination of time-dependent hat functions P_{k}(r,t)
(similar to the approximation with static hat functions in
Eq. 6)

(11)

Here, the hat functions are - like in the Claverie approach - step-wise defined triangles

(12)

but with grid points r_{k }that change with time as

(13)

where the position r_{k}(t) is that of an ideal sedimenting particle
with an s-value s_{G} . The exponential movement of the
grid-points is abbreviated as a(t-t_{0}).
Only the meniscus *m* and bottom *b* are fix, since they have to denote
the constant solution column position measured as distance from the center of
rotation.

There is one obvious problem - how can we choose a good set of radial grid points for the entire simulation, when they permanently 'sediment' towards the bottom of the simulated centrifuge cell? The key is to use a properly spaced initial grid consistent with the movement of grid points, and then to change the grid at well-defined times. This works as follows:

If we choose initially the radial grid exponentially as

(14)

combined with the movement in time described by Eq. 13, we find that after a well-defined time

(15)

the grid is mapped precisely onto itself:

(16)

This is a very useful property. Because (except for the meniscus and bottom) every grid point is now exactly where it's right-hand neighbor was before, we change the grid by removing the point N-1, and inserting a new second point (designating the previous second point now as #3). This is illustrated in the following cartoon:

In this way, if our time-increments are exactly Dt_{swap},
all grid points remain virtually stationary except for their index!.
Therefore, the movement can be accomplished
by simply *renumbering* all indices from 2 to N-2. (This procedure is
similar to the phenomenon when looking at a spinning wheel with a stroboscopic
lamp that has the same frequency as the rotation.)

As described in (Schuck, 1998), the complete algorithm can be mathematically
put together in a form very similar to the Claverie algorithm described
above: Multiplication of the Lamm equation with our (now time-dependent)
hat function P_{k}(r,t) and partial integration leads to

(17)

Insertion of our approximation of c(r,t) as superposition of time-dependent hat functions (Eq. 11, with the running index now called 'j') gives us expressions for the matrix elements

(18)

New here (in comparison with Eq. 8a) is the *time-derivative* of P_{j}(r,t),
which can be easily calculated as

(19)

(except for the first and last hat function, which has one leg that stays
constant in time). At this point, we can transform Eq. 18 by change of variable r(r,t)
= r/a(t-t_{0}), which corresponds to the transition to
the moving frame of reference. This allows us to separate the
time-dependence of our hat functions P_{k}(r,t) out of the
integrals. Analogous to the case of the Claverie approach, the resulting
integrals can be calculated analytically and are abbreviated as **A**^{(1)},**A**^{(2)}
, **A**^{(3) }and **B**. Their numerical values can be found
here (in C-code). Since we have the constant position of meniscus and bottom, the corner 2x2
sub-matrices are time-dependent, but all other elements are not.

Using the matrix notation, we can simplify Eq. 18 to

(20)

which is very similar to Eq. 8b in the Claverie approach. There are two
important differences: 1) there is the factor a(t-t_{0}) which takes into
account the changing diffusion distances in the moving reference frame (because
of the stretching of the reference frame); 2) there is an extra term with **A**^{(3)},
which was introduced by the time-derivative of the hat-functions.
Importantly, it can be shown that **A**^{(2)}
- **A**^{(3) }= 2**B**, which *eliminates the sedimentation term
for s = s _{G}*, i.e. when the frame of reference sediments with the
same s-value as the particle. The remaining expression from combination of

In summary, we have formulated a finite element solution that essentially
requires only diffusion and radial dilution to be calculated if we use *s = s _{G}*.
As a
consequence, the limiting case of no or small diffusion is trivial to solve with
this algorithm. On the other hand, if we choose

Similar to the other approaches described above, Eq. 20 can be used with higher precision in the Crank-Nicholson scheme. For details see Schuck (1998).

**The implementation in Sedfit **

Because the moving grid simulation requires time-steps not larger than Dt_{swap
}given by Eq. 15, the moving grid with *s = s _{G}* is only
used for cases with

Several control parameters for the simulation can be changed in Sedfit (e.g., the number of radial grid
points, a tolerance parameter for the adaptive step-size; see the parameter box).
Also, for the interested user, the different algorithms can be chosen with
different parameters, fixed Dt, *s _{G}*,
etc, in order to enable exploration of the effect of these parameters on the
execution time and precision.

In both algorithms, changes in the rotor speed with time can be simulated by simply calculating the Eqs. 20 and 9b using a different average rotor-speed for each time-interval Dt. This is useful for simulating the acceleration phase of the experiment.

As a variant, a linearized version of the Lamm equation can be solved, which is useful when working with external loading cells. Further, a linear Lamm equation with constant force (as an electrophoresis model) has been implemented. In both cases, the algorithm is the same, only requiring exchange of the matrix elements.

Rapid self-association of the sedimenting macromolecule (i.e. self-association where the local re-equilibration of the oligomeric species is fast compared to the sedimentation and diffusion process) is simulated by using locally weight-average sedimentation coefficients and gradient average diffusion coefficients, as described first by Cox (1965), but implemented into the finite element algorithms. This requires additional development of s and D into hat-functions:

If these sums are inserted and carried through the same finite element formalism, all the integrals will contain three hat functions with three different indices, and the matrix equations expand to tensor equations. These can still be calculated very efficiently, in particular since the Crank-Nicholson scheme can be replaced by a more general predictor-corrector scheme (Schuck, 1998).

Similarly, non-ideality in form of concentration-dependent s(c) and D(c) is implemented analogous to the discretization in Eq. 21.

The extension of the algorithm to heterogeneous associations (A + B <-> C) requires the solution of a set of three Lamm equations, coupled with chemical reaction terms. The finite element solution of the coupled system has been derived, and is currently in the implementation stage.

Note, as of 2007,a significant update and improvement in the Lamm equation solution is implemented, as described in

**P.H. Brown, P. Schuck (2007) A new adaptive
grid-size algorithm for the simulation of sedimentation velocity profiles in
analytical ultracentrifugation. **

**Click on this link to read how this works: **
Computer Physics Communications in
press

**This work is not yet included in this introductory page.**

**References:**

Behlke,
J., and O. Ristau. 1997. Molecular Mass Determination by Sedimentation Velocity
Experiments and Direct Fitting of the Concentration Profiles.* Biophys. J. *72:428-434

Bethune,
J.L., and G. Kegeles. 1961. Countercurrent distribution of chemically reacting
systems. I. Polymerization. *J. Phys. Chem.* 65:433-438

P.H. Brown, P. Schuck
(2007) A new adaptive grid-size algorithm for the simulation of sedimentation
velocity profiles in analytical ultracentrifugation.
Computer Physics Communications in
press

P.H. Brown, P. Schuck (2006) Macromolecular size-and-shape distributions by
sedimentation velocity analytical ultracentrifugation.
Biophysical Journal 90:4651-4661

Cann,
J.R., and G. Kegeles. 1974. Theory of Sedimentation for kinetically controlled
dimerization reactions. *Biochemistry* 13:1868-1874

Claverie,
J.-M., H. Dreux, and R. Cohen. 1975. Sedimentation of generalized systems of
interacting particles. I. Solution of systems of complete Lamm equations. *Biopolymers.*
14:1685-1700.

Cox,
D. J. 1965. Computer simulation of sedimentation in the ultracentrifuge. II.
Concentration-independent sedimentation. *Arch. Biochem. Biophys. *112:259-266.

J. Crank and P. Nicholson. 1947 A practical method for numerical evaluation of solutions of partial differential equations of the heat-conducting type. Proc. Cambridge Philos. Soc. 43:50-67

J. Dam, C.A.
Velikovsky, R.A. Mariuzza, C. Urbanke, P. Schuck (2005) Sedimentation velocity
analysis of heterogeneous protein-protein interactions: Lamm equation modeling
and sedimentation coefficient distributions c(s).
__
Biophysical Journal__
89:619-634

Dishon,
M., Weiss, G.H., and D.A. Yphantis. 1966. Numerical solutions of the Lamm
equation. I. Numerical procedure. *Biopolymers* 4:449-455.

Durran, D.R. Numerical methods for wave equations in geophysical fluid dynamics. Springer, New York, 1999, p. 139

Faxén,
H. 1929. Über eine Differentialgleichung aus der physikalischen Chemie. *Ark.
Mat. Astr. Fys. *21B:1-6.

Fujita, H. 1962. Mathematical Theory of Sedimentation Analysis. Academic Press, New York

Holladay,
L.A. 1979. An approximate solution of the Lamm equation. *Biophys. Chem. *10:187-190Lamm,
O. 1929. Die Differentialgleichung der Ultrazentrifugierung. *Ark. Mat. Astr.
Fys. *21B:1-2

Marque,
J. 1992. Simulation of the time course of macromolecular separations in an
ultracentrifuge. II. Controlling the solute concentration. *Biophys. Chem. *42:23-27

Minton,
A.P. 1992. Simulation of the time course of macromolecular separations in an
ultracentrifuge. I. Formation of a cesium chloride density gradient at 25ºC. *Biophys.
Chem. *42:13-21

Philo,
J.S. 1997. An improved Function for Fitting Sedimentation Velocity Data for
Low-Molecular-Weight Solutes. *Biophys. J. *72:435-444

P. Schuck,
C.E. McPhee, and G.J. Howlett. (1998) Determination of
sedimentation coefficients for small peptides. * Biophysical Journal*
74:466-474

P. Schuck (1998) Sedimentation analysis of non-interacting and
self-associating solutes using numerical solutions to the Lamm equation. * Biophysical
Journal*
75:1503-1512