- Email: [email protected]
Computer simulation of sedimentation in the ultracentrifuge
ARCHIVES
OF
BIOCHEMISTRY
Computer
AND
BIOPHYSICS
Simulation
112, 249-258 (1965)
of Sedimentation
in the
Ultracentrifuge
I. Diffusion DAVID Clayton
Foundation
Biochemical
J. COX
Institute and the Department of Texas, Austin, Texas Received
of Chemistry,
The University
June 4, 1965
The procedure of Vink for the computer simulation of concentration-independent diffusion in a rectangular cell is evaluated. The results of simulated experiments are found to be in excellent agreement with those derived from the exact solution of Fick’s second law. The Vink model is modified to provide methods for the simulation of concentration-independent diffusion in a sector-shaped cell and of concentrationdependent diffusion. The results of representative calculations are presented and discussed.
The exact solution of the continuity equation describing the sedimentation in the ultracentrifuge of any but the most simple systems presents very great mathematical difficulties (1). It may be of interest to predict the sedimentation behavior of systems for which analytic solutions are not available. Such predictions might be possible if procedures could be developed for the simulation of the sedimentation process by a sequence of discrete steps described in terms suitable for the application of a digital computer. The ingenious counter-current distribution analog of Bethune and Kegeles (2) has proved useful for the prediction of the shapes of the sedimenting boundaries of reacting systems. This important method is readily applied to the sedimentation, in a rectangular cell under a constant field, of components whose sedimentation and diffusion coefficients are independent of concentration; the extension of the method to the concentration-dependent case in a sectorshaped cell in a field which varies with position seems quite difficult. Moreover, the effect of the use of a discontinuous model on the accuracy of the results is not readily evaluated. For these reasons, it may be fruitful to explore other simulation procedures more closely analogous to the actual
physical processes occurring in a sedimenting system. In this and the following paper, a procedure is developed which considers sedimentation and diffusion separately. The components of the system are said to sediment without diffusing for a short time (At) and then to diffuse without sedimenting for an equal time. If the time interval is short, the procedure should adequately simulate real systems in which sedimentation and diffusion occur simultaneously. The advantage of dealing with sedimentation and diffusion separately is that, as will be seen, the inadequacies of a proposed computational method can be more easily detected and eliminated. If the procedure is to be useful, a large number (k) of successive alternate rounds of simulated sedimentation and diffusion should reproduce the actual physical situation at a time k At after the establishment of the specified initial conditions. The effect of the use of a discontinuous model may be tested in two ways. First, a relatively simple case may be selected for which an analytic solution is available. A concentration distribution or boundary shape is calculated using the discontinuous model to be tested. The result is then compared with the exact solution. For the more
249
250
DAVID
interesting cases where the correct result is not known, several separate model calculations can be done, dividing the total time of the experiment to be simulated into different numbers of time intervals; the impact of the discontinuous nature of the description can be evaluated by determining to what degree the results are altered by the arbitrary selection of the number of steps into which the total process is divided. During an ultracentrifuge run, a sedimenting boundary spreads as a consequence of the diffusion of the sedimenting component. In order to simulate such an experiment, it is therefore necessary to develop a computational model for diffusion in the sector-shaped cells generally used with the ultracentrifuge. Vink (3) has reported a simulation procedure for concentrationindependent free diffusion in one dimension in a rectangular cell. In the present paper, the accuracy of Vink’s model is evaluated by comparing the results obtained by its use with the exact solution of Fick’s second law. The method is then suitably modified to deal with cases where diffusion occurs in a sector-shaped cell or where the diffusion coefficient varies with concentration.
J. COX
continuous
analog of Fick’s first law:
SIMULATION AND
PROCEDURES RESULTS
Concentration-independent diffusion in a rectangular cell. The simulation procedure for this simplest case has been published by Vink (3). The description of the method is repeated in detail here to make intelligible the modified procedures described later in the paper. The rectangular cell, with crosssection A, is divided into n compartments by n + 1 boundaries equally spaced along the long axis of the cell; thus, box i is enclosed by boundaries i and i + 1. The interval between the boundaries is Ax, and the distance between the midpoints of adjacent boxes is also Ax. The flow of mass into box i from box i - 1 (AmI) during a short time interval At is given by a dis-
Ci) )
(la)
where D is the diffusion coefficient and Ci is the concentration in compartment i. The mass flow into box i from box i + 1 (Amz) is given by the similar expression:
Am;= D.A.&
(“+’ - “) . (lb) Ax
The sum of Am, and Am, is the change, during At, in the mass of the diffusing component in box i. Normally, one of the Am terms is positive and the other is negative; the sum of Am, and Am2 may be either posi-tive or negative. The change in the concentration of the diffusing component in box i is: ACi
=
AmI + Amz = D.At. A.Ax
;
(&+I - Ci) _ (Ci - ci-1) Ax
-
Ax:
1
@a)
Equation (2a) will be recognized as a discontinuous “translation” of Fick’s second law :
dc -=Dol”c
METHODS Programs incorporating the simulation procedures to be described were written in FORTRAN63 and the calculations were done by a Control Data Corporation 1604 digital computer.
a-‘,,
= D.A.At
Anal
dt
Equation putational ACi
(2b)
dx2 *
(2a) can be simplified purposes: =
,-~(Ci+l
+
ci-1
-
for com-
2ci)
D-At ff=(az)z* Since there can be no mass flow across boundaries 1 and n + 1 at the ends of the cell, the appropriate expressions for compartments 1 and n are: AC1 = a(& - C,)
(24
AC,
(24
= LY(&~
-
c,).
A diffusion experiment is simulated by the following procedure. Given the diffusion coefficient, a suitable box dimension and time interval are selected, and a! is calculated. The initial concentration distribution is specified; for example, if there is initially a
SIMULATION
sharp boundary between solvent and solution at the center of a lOO-box system, the distribution is described by: c; = 0, i = 1;**50 (3)
ci = co, i = 51;.*100,
where Co is the initial concentrat.ion on the solution side of the boundary. The procedure is not restricted to t.his one initial boundary shape. With the initial concentrat.ions specified, Eqs. (2c, 2d, and 2e) arc applied to each compartment. When the calculation has been done for each box, a new set of conccntrations ir calculated. Ci’ = Ci + ACi.
(4)
The set of C’i values defines the concentration distribut.ion after time At. The new concent,rat,ions arc insert.ed into Eqs. (2c, 2d, and 2e), and the calculation is repeated to give the distribution after time 2eAt. Afler h: such transfers, the calculated concentration distribution should reproduce that which would have been observed in a real experiment. k At aft.er the establishment of the initial conditions. The results may bc conveniently converted to a concent,ration gradient, distribution :
0
dc (Ci - Ci-1) d; i = Ax
’
251
OF DIFFUSION
(5)
The gradient subscript in Eq. 5 refers to boundary i; the concentration subscript indicates box i. Gradients are defined only for boundaries 2 through n. Vink (3) has reported using the above procedure for several t.ransfers with a single value of cr; the results appeared promising, in that the calculated boundaries were symand at. Icast approximately metrical Gaussian. The practical use and further elaborat,ion of the method require, however, a detailed examination of the accuracy of the results. In particular, it is essential to know whether the choice of (Y is completely unrestricted. It must be known whether, with a given D, the box interval Ax and the time per t.ransfer At can bc assigned any arbitrary values. It can be shown immediately that. cy is
restricted to values below 0.5. This restricbion is readily understood by considering what happens during the ‘first. t.ransfer in a system with a sharp initial boundary between solvent in boxes 1 through i and solution in boxes i + 1 through n. Examination of Eq. (2~) reveals t.hat, after the first transfer with (Y greater than 0.5, the concentration in box i will be greater t.han that in box i + 1. That is, the concent.ration gradient at the original boundary will be reversed in sign, and a peak in the conccntrat,ion distribution will result. If CYis greater than 1.0, then, after the first transfer, the calculated concentration in box i + 1 will be negative. Thus, with a given diffusion coefficient, physically real results certainly cannot be obtained unless Ax and At arc so chosen that DAM/* is less than 0.5. It remains to see how much effect variation of (Y within the range below 0.5 has on the calculated distribution. A given system may be divided into a large number of narrow compart,ments or a small number of wide boxes; or a given total time Ii A1 may be described by a large k and a small At or vice versa. The fundamental question is whct.her the physical process is accurately depicted with any values of a! and hi. This question was approached empirically by applying the procedure to the diffusion of a component wit.h D = 8.0 X lo-’ cm* per second in a system of 100 boxes, each 0.01 cm wide. At zero time there was a sharp boundary between the solvent and a solution whose concentration was 1.0 mg per milliliter. The gradient curve after 5000 seconds was calculated for this system. Since &
= k $f2
= ka = 40.0,
(6)
the simulation calculation was done by using four different pairs of k and (Y, each selected so that the product of the two was 40.0. In addition, t.hc correct shape of the gradient was calculated by using the exact solut.ion of the appropriate continuity cquation : dc -= dx
co exp (-x2/4Dt). 2(7rDL)“2
(7)
252
DAVID TABLE
SIMULATED
DIFFUSION
I
IN A RECTANGULAR
Diffusion parameter:
0.025
0.05
0.10
0.20
Nzber transfers:
1600
800
400
200
11%- %I1 (cm) 0.00
0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20
CELL
Exact Solution
Concentration gradient (mg/cm~)
4.47 4.36 4.04 3.56 2.99 2.38 1.81 1.31 0.90 0.59 0.37
4.47 4.35 4.04 3.56 2.99 2.38 1.81 1.31 0.90 0.59 0.37
4.46 4.35 4.04 3.56 2.99 2.39 1.81 1.31 0.90 0.59 0.37
4.46 4.35 4.04 3.56 2.99 2.39 1.81 1.31 0.90 0.59 0.37
4.46 4.35 4.04 3.56 2.99 2.39 1.81 1.31 0.90 0.59 0.37
The calculated gradient distributions are displayed in Table I; the results are clear and somewhat surprising. Variation of (Y over a eightfold range from 0.025 to 0.20 has essentially no effect on the shape of the gradient curves. Moreover, the results of the simulation technique are in excellent agreement with the exact solution for all values of ar; the small differences observed are within the range of errors normally encountered in reading t,he schlieren record of a diffusion experiment. The conclusion may be drawn that any convenient choice of Ax and At will be acceptable for the simulation of diffusion by the Vink procedure if the resulting value of (Y is between 0.025 and 0.2; the choice may in fact be more flexible than this. Variation of (Y within me calculation. It is now necessary to modify the simulation procedure described above in order to make it suitable for describing the two cases of immediate interest : concentration-dependent diffusion and diffusion in a sectorshaped cell. Although the specific procedures to be used in these two cases are quite different, the two have one feature in common: the diffusion parameter a! must vary from point to point in the cell. In the sector cell, the cross sections of the successive boundaries and the volumes of the successive boxes increase outward from the axis of the sector. When the diffusion coefficient depends on the local concentration of the solu-
J. COX
tion, and the concentration changes across the cell, then D and thus (Y must vary from point to point. Before proceeding to deal with either case, it is necessary to find out whether the use of different diffusion parameters in different parts of the cell will, of itself, deform the calculated boundary. The simplest way to test this possibility is to recast the rectangular-cell, constant-D case in such a way that (Yvaries from one compartment to the next; this can be done by dividing the cell into compartments of different instead of equal widths. A series of n + 1 boundary positions is specified; an expression analogous to Eq. (2a) is derived for the change, during the time interval At, of the concentration in box i, enclosed between boundaries i and i + 1: DAt Aci = (Xi+1 - Xi)
(8) (Zi - f&l) (is+1 - 5%) 1.
cc,+1 - CJ [
(Ci - ci-1)
In Eq. (8), xi is the position of boundary i; Zi is the position of the midpoint of box i. It should be stressed that the assignment of different widths to the compartments is entirely artificial; the physical process being simulated is still concentration-independent diffusion in a rectangular cell. The curve plotted from the results should therefore be identical to one obtained using a constant box interval; in particular, the calculated gradient curve should be symmetrical and Gaussian. Table II shows the result of one such test. A lOO-box cell was constructed as follows: the width of box 1 was taken as 0.01 cm. The widths of the other boxes were specified by:
(xi+1 - Xi) = (Xi - Xi-l) . g.
(9)
Thus, the width of box i was 0.01 (6.01/ 6.00)i-1, and the boxes increased progressively in width from 0.0100 to 0.0118 cm. Since (Y is inversely proportional to (Ax>~, a! at box 1 was about 1.4 times its value at box 100. The simulation procedure of Eq. 8 was applied for 1000 transfers with DAt equal to 5 X lo+. The resulting gradient distribution was tested for sym-
SIMULATION TABLE
II
TABLE
DIFFUSION SIMULATION WITH CONTINUOUS VARIATION OF Box DIMENSION, TEST FOR SYMMETRY OF GRADIENT CURVE” Box boundaries in lower limb lb - all (4
Concentration gradients in lower limb (mg/cm4) Observed
0.00000 0.01085 0.02170
3.994 3.970 3.900
0.03260 0.04355 0.05450 0.06545 0.07640 0.08735 0.09835
3.786 3.631 3.440 3.221 2.979 2.722 2.458
0.10940 0.12045 0.13150 0.14260 0.15370 0.16480 0.17595 0.18710 0.19825 0.20945 0.22065
2.191 1.929 1.679 1.442 1.224 1.026 0.849 0.694 0.560 0.446 0.350 0.272 0.208 0.158 0.117 0.087 0.063 0.045 0.032 0.023 0.016
OF DIFFUSION
Calculated by interpolation in upper limb
3.994 3.970 3.900 3.786 3.631 3.441 3.221 2.980 2.722 2.457 2.190 1.929 1.677 1.441 1.222 1.024 0.847 0.692 0.558 0.444 0.349 0.271 0.208 0.157 0.117 0.087 0.063 0.045 0.032 0.023 0.016
DISCONTINUOUS Box number 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
III
CEANQE IN Box
DIMENSION
Concentration @g/ml) Variable AZ
Constant Ax
1.000 0.999 0.999 0.998 0.997 0.996 0.993 0.989 0.984 0.977 0.966 0.953 0.935
1.000 0.999 0.999 0.998 0.997 0.995 0.993 0.989 0.984 0.976 0.966 0.952 0.934
0.912
0.911
0.8&l 0.850 0.809 0.763 0.711 0.655 0.594 0.532
0.883 0.849 0.808 0.762 0.711 0.654 0.594 0.532
this case, the box dimension was 0.01 cm for all compartments on the solution side of the 0.24315 initial boundary and 0.02 cm for all compart0.25440 ments on the solvent side. There was a 0.26570 sudden doubling of the box dimension and, 0.27700 thus, a fourfold change in (Y at the center of 0.28830 the diffusing boundary. The concentration 0.29965 0.31100 distribution on the solution side of the 0.32235 boundary after 200 transfers with D At 0.33375 = 1.0 X 10m6is shown in Table III; the initial concentration was 1.0 mg per millidifa Simulated concentration-independent liter. The distribution calculated with the fusion in a rectangular cell. Compartment disame parameters in a set of compartments mensions selected as described in the text. 1000 0.01 cm wide on both sides of the boundary transfers with DAt = 5 X 10-B. is shown for comparison. It is apparent t’hat the correct result is obtained for the recmetry by reflecting the upper limb of the curve through the original boundary posi- tangular-cell, constant-D case, even when (Y is arbitrarily assigned different values in diftion and examining the result to see if it ferent parts of the cell. Simulation procecoincided with the lower limb of the observed curve. As Table II indicates, the dures may now be developed for diffusion in gradient curve was almost perfectly sym- a sector cell and for concentration-dependent diffusion. The results of Tables II and metrical. The sensitivity of the rectangularcell, constant-D simulation to variation in cr III lend some confidence that the distribuwas subjected to a more rigorous test by de- tions calculated for these cases will not be distorted by the fact that the diffusion pafining still another set of compartments. In 0.23190
254
DAVID
rameter happens to vary with position in the cell. Diflusion in a sector-shaped cell. The physical situation to be described is the spreading of a boundary in a synthetic boundary cell in an ultracentrifuge operated at low speed. In the usual notation, the sector angle is 8, and the height of the sector is b. The cell is divided into n compartments by n + 1 cylindrical surfaces located at various distances r from the axis of the sector; boundary 1 is the meniscus, and boundary n + 1 is the cell bottom. The radial dimensions of the various compartments need not be the same. The area of boundary i is :
(10)
Ai = mi.
The volume of box i, between boundaries i and i + 1, is: vi = W&i+1
- ri),
(11)
where pi is the radial position of the midpoint of box i: i;,
=
b-i
z
+
ri+11
2
(12)
*
During the time At, the flow of mass into box i from box i - 1 across boundary i is:
and the flow of the diffusing component into box i from box i + 1 across boundary i+ 1 is: (Cc+1 - Cd Am2 = D. At. bf?ri+l (pi+l _ pi) .
(13b)
The change in the concentration of the diffusing component in box i during time interval At is: AC,
z-
_
Am1 + Am2 Vi
=
J. COX TABLE IV COMPARISON OF SIMULATED DIFFUSION RECTANGULAR AND SECTOR-SHAPED CELLS
Concentration gradient (mgjcm’)
Since the box dimensions during the calculation, the putations may be reduced diffusion parameters bi and
1
(14)
do not change number of comby defining the fi:
Rectangular cell
sector cell
0.508 0.625 0.762 0.918 1.094 1.290 1.504 1.735 1.980 2.234 2.493 2.750 3.000 3.235 3.449 3.636 3.788 3.901 3.970 3.994 3.970 3.900 3.786 3.631 3.441 3.221 2.980 2.722 2.457 2.190 1.929 1.677 1.441 1.222 1.023 0.847 0.692 0.558 0.444
0.516 0.634 0.773 0.930 1.107 1.305 1.521 1.753 1.998 2.253 2.511 2.768 3.017 3.251 3.464 3.647 3.797 3.907 3.973 3.993 3.967 3.893 3.776 3.619 3.427 3.206 2.962 2.704 2.439 2.172 1.911 1.660 1.425 1.207 1.011 0.835 0.682 0.549 0.438
6.3126 6.3231 6.3337 6.3442 6.3548 6.3654 6.3760 6.3866 6.3973 6.4079 6.4186 6.4293 6.4400 6.4508 6.4615 6.4723 6.4831 6.4939 6.5047 6.5155 6.5264 6.5373 6.5482 6.5591 6.5700 6.5810 6.5919 6.6029 6.6139 6.6250 6.6360 6.6471 6.6581 6.6692 6.6804 6.6915 6.7026 6.7138 6.7250
DAt c;i(ri+1 - f-i>
_ ri+1cci+1- Ci> (Pg.1 - Pi)
IN
bi = I
Ji =
DAt ri+l Pi(Ti+1 - ri)(fi+l _ Ti(Ti+l
thus simplifying ACi = fi(C,,
- Pi) (15)
DAt r; -
ri)(Pi
-
f&l)
’
Eq. (14) to give: - Ct> + b;(Ci+l -
Ci). (16)
SIMULATION
The condition that no mass flows across boundaries 1 and n + 1 at the ends of the cell is: fl = b,+1 = 0.
255
OF DIFFUSION
(17)
Diffusion in a sector-shaped cell is simulated in the following way. A convenient set of boundary locations ri is specified. The locations of the midpoints of the boxes, pi, are calculated. A time interva.1 At is selected and, with D given, the diffusion parameters bi and fi are determined for each box; At must be so chosen that each of the diffusion parameters is less than 0.5, and preferably lies within the range 0.025-0.2. The initial concentration distribution is specified, and Eq. (16) is applied to each box. The concentration distribution at time At is calculated, and the process is repeated k times. The calculated concentration distribution should reproduce the physical situation k At after the beginning of a real experiment. Table IV shows gradient distribution resulting from such a calculation. The set of boundaries chosen for the simulation happened to be the ones specified by Eq. (9) ; i.e., the width of the boxes varied across the cell. Thus the results can be compared directly with those shown in Table II. There is, of course, no necessity for specifying a variable cell dimension in order to describe the sector-cell case. The model cell was placed between 6.00 and 7.08 cm from the axis of the sector; the model was thus very much like the usual ultracentrifuge cell. As in the experiment of Table II, the initial sharp boundary was placed between the 50th and 51st boxes; the initial concentration on the solution side of the boundry was; 1.0 mg per milliliter, and diffusion was simulated by 1000 transfers with D At equal to 5 X 10-6. Inspection of Table IV reveals that the effect ot the sector shape of the cell is small but noticeable. The shape ot the gradient curve is not detectably changed, but the whole curve is shifted very slightly (about 7~) toward the narrow end of the sector. Concentr&on-dependent diffusion in a rectangular cell. In simulating this important case, the diffusion coe@icient at infinite dilution is specsed, and the dependence of D
on concentration series :
is described by a power
D = DO (1 + k# + kzC2 +. . . -).
(18)
The concentration which determines the local diffusion coefficient at boundary i is taken to be the mean of the concentrations on the two sides of the boundary, in boxes i and i - 1. Diffusion where D is a linear function of concentration has been considered in detail; in this case, the local diffusion coefficient at boundary i is: Di = Do 1 + k (ci +2 “-‘))
.
(19)
It is convenient to divide the cell into n compartments by n + 1 equally spaced boundaries, as in the concentration-independent case. Equation (2a) is modified to yield a working equation suitable for concentration-dependent diffusion: 1+
ACi = 00
k(Ci
+
Ci-1)
>
K
* (Ci-1 - C;) +
1 + ,“(” * (&+I
: “‘+l)> -
Ci>
(2Oa)
1 7
where
(20b) The fact that no mass crosses the ends of the cell, boundaries 1 and n + 1, is stated as follows : k(Cl + 2
c2,>
(2Oc)
* cc2 - Cd G’Od)
*(C?r,- cn>. A concentration-dependent diffusion experiment was simulated with the working Eqs. (20a), (~OC), and (20d); the result is shown in Fig. 1. The model cell consisted of 100 boxes 0.01 cm wide. The sharp boundary
256
DAVID
J. COX
50.0 7 E 9 40.0 If ‘: .? 30.0 ‘0 E!
W
.z 20.0 ‘0 L ‘; 0” 6 10.0 0
-3.0
-2.5
-2.0
-1.5
-1.0
-0.5
0
X-X,
to.5
t1.0
t1.5
t2.0
+2.5
+3.0
3.5
(mm)
FIG. 1. Simulation of diffusion with D a linear function of concentration. Diffusion coefficient varies from 2.5 X lO+ cm2/second at C = 0 to 7.5 X 10-7 cm2/second at C = 10 mg/ml. Triangles show the concentration gradient distribution after 8000 seconds. Circles show the gradient distribution resulting from the diffusion of a component with D = 5.0 X 10-7 cmz/second for 8000 seconds.
between the solvent and a solution at 10 mg per milliliter was formed between boxes 50 and 51. For illustrative purposes, an extreme case was selected in which the diffusion coefficient varied linearly with concentration from 2.5 X lo-’ cm2 per second at infinite dilution to 7.5 X lO-’ cm2 per second at 10 mg per milliliter. The time interval per transfer was 10 seconds, and the simulation was carried out for 800 transfers. The diffusion for 8000 seconds of a component with D equal to 5.0 X lo-’ cm2 per second at all concentrations was also simulated; the resulting gradient curve is shown in Fig. 1 for comparison with the analogous concentration-dependent experiment. Gosting and Fujita (4) have reported an exact solut’ion of the continuity equation for the case where the diffusion coefficient is a function of concentration. The expression obtained for the concentration gradient distribution contains a number of terms in the form of a power series in initial solute concentration. Costing (5) has shown a plot based on this solution which is analogous to Fig. 1; it compares the gradient curves for concentration-dependent and concentrationindependent diffusion. However, the plot for the concentration-dependent case is an approximation; it takes account of only the
TABLE COMPUTER
V
TIME REQUIRED FOR VARIOUS SIMULATION PROCEDURES
Process simulated
NO. transfers
Rectangular-cell, constant D Sector-cell, constant D Rectangular cell D - Do(1 + kc) Rectangular cell D - Do(1 + MY)
1000
1 minute
14 seconds
1000
1 minute
20 seconds
1000
1 minute
35 seconds
500
3 minutes
17 seconds
Computer time
first two terms of the exact solution. Fig. 1 and Gosting’s Fig. 6 are quite similar. The gradient distribution for the concentrationdependent system is markedly skewed toward the region of the lower diffusion coefficients. The plot for the concentrationdependent case crosses the analogous curve for the system with a constant diffusion coefficient at three points. In the graph presented by Gosting, the crossover points occur at the peak of the symmetrical concentration-independent plot and at two other points equidistant from the maximum of the constant-D gradient distribution. In Fig. 1, the middle crossover point is dis-
SIMULATION
placed from the peak of the Gaussian constant-D plot toward the region of lower concentrations, and the other two crossover points are not symmetrically arranged around the peak. The differences are probably due to the neglect of the higher terms of the analytic solution in the construction of Gosting’s Fig. 6. The sum of these higher terms is apparently negative and approxirnately symmetrical (4). No great amount of computer time is needed for calculations of this kind. Table V shows the time required to process simulation programs of several types. DISCUSSION
The procedures described in this paper were devised as a part of a more general project leading to the development of methods for the simulation of sedimentation in the ultracentrifuge. In this context, the major results of the work are that the Vink procedure describes diffusion very accurately, and that the method can easily be modified to deal with diffusion in a typical ultracentrifuge cell. The simulation of diffusion may be of some practical use in itself. Exact solutions of Fick’s second law are, of course, available and are widely used in systems where the variation of the diffusion coefficient with concentration is negligible; the simulation procedure has no important practical use in these simple cases. Where the diffusion coefficient varies with concentration, however, the methods presently in use leave something to be desired. Exact solutions are available for handling the results of diffusion experiments when D is a function of solute concentration (4). These exact solutions are somewhat cumbersome in practice, however; their use involves a great deal of laborious computation. As a consequence, methods for handling such data generally rely on approximations to the exact solution. The most commonly used simplification consists of discarding the higher terms in the series solution (6). The total effect of such a mathematical approximation on the accuracy of the calculated results is generally not obvious, but one assumption implicit in this procedure is that the dependence of D
OF DIFFUSION
257
on powers of concentration above the first is negligible. When diffusion coefficients at different solute concentrations are calculated using the standard approximate treatments, it is not unusual to find that D does not vary linearly with concentration. Infinite-dilution diffusion coefficients derived by extrapolation from such data are of uncertain accuracy, particularly if the concentration dependence of the diffusion coefficients happens to be steep. The simulation technique may be suggested as a method which would place the calculation of diffusion coefficients for concentration-dependent systems on firmer ground. An infinite-dilution diffusion coefficient and the parameters describing the dependence of D on concentration might be calculated by the usual approximate methods. A diffusion experiment could then be simulated by using these approximate parameters and the result of the simulation could be compared with the observed boundary shape. If the calculated and observed boundaries did not coincide, the diffusion parameters could be adjusted to give a better fit. In principle, there is no reason why the computer could not be presented with the experimental concent,ration or gradient distribution and instructed to find a Do and regression coefficients that would generate the observed results. However, a program of this type has not yet been written. The most attractive feature of the simulation technique is its great flexibility. The exact solution of the continuity equation is very much complicated by a two- or threeterm concentration dependence of the diffusion coefficient, but little more computer time is required for the description of such a system by the simulation technique t,han for the simple constant-D calculation. Moreover, the initial boundary in a diffusion experiment is usually not perfectly sharp. The methods generally used for diffusion calculations formally handle a blurred initial boundary in terms of a time error. This procedure contains an implicit assumption about the shape of the initial boundary; for example, it assumes, for the constant-D case, that the initial boundary is Gaussian.
258
DAVID J. COX
The application of t.he simulation technique requires only that. the shape of t.he initial boundary be known, so that it may be specified at the beginning of the calculation. REFERENCES Theory of Sedi1. FUJITA, II., “Mathematical mentation Analysis.” New York, Academic Press (1962).
2.
J. L., AND KEGELES, G., J. Phjys. Chem. 66, 1761 (1961). 3. \'Ih'K, II., nCk& Chem. &and. 18, 409 (1964). 4. GOSTIN, L. J., ASD FLTJIT.~,I-I., J. Am. Chem. sot. 79.1359 (1957). 5. (doSTISG, 1~. J., A&an. Protein Chem. 11, 429 (1956). 6. CHEETH, J. M., J. Am. Chem. Sot. 77. 64% (1955). BE?.IIUNE,
OF
BIOCHEMISTRY
Computer
AND
BIOPHYSICS
Simulation
112, 249-258 (1965)
of Sedimentation
in the
Ultracentrifuge
I. Diffusion DAVID Clayton
Foundation
Biochemical
J. COX
Institute and the Department of Texas, Austin, Texas Received
of Chemistry,
The University
June 4, 1965
The procedure of Vink for the computer simulation of concentration-independent diffusion in a rectangular cell is evaluated. The results of simulated experiments are found to be in excellent agreement with those derived from the exact solution of Fick’s second law. The Vink model is modified to provide methods for the simulation of concentration-independent diffusion in a sector-shaped cell and of concentrationdependent diffusion. The results of representative calculations are presented and discussed.
The exact solution of the continuity equation describing the sedimentation in the ultracentrifuge of any but the most simple systems presents very great mathematical difficulties (1). It may be of interest to predict the sedimentation behavior of systems for which analytic solutions are not available. Such predictions might be possible if procedures could be developed for the simulation of the sedimentation process by a sequence of discrete steps described in terms suitable for the application of a digital computer. The ingenious counter-current distribution analog of Bethune and Kegeles (2) has proved useful for the prediction of the shapes of the sedimenting boundaries of reacting systems. This important method is readily applied to the sedimentation, in a rectangular cell under a constant field, of components whose sedimentation and diffusion coefficients are independent of concentration; the extension of the method to the concentration-dependent case in a sectorshaped cell in a field which varies with position seems quite difficult. Moreover, the effect of the use of a discontinuous model on the accuracy of the results is not readily evaluated. For these reasons, it may be fruitful to explore other simulation procedures more closely analogous to the actual
physical processes occurring in a sedimenting system. In this and the following paper, a procedure is developed which considers sedimentation and diffusion separately. The components of the system are said to sediment without diffusing for a short time (At) and then to diffuse without sedimenting for an equal time. If the time interval is short, the procedure should adequately simulate real systems in which sedimentation and diffusion occur simultaneously. The advantage of dealing with sedimentation and diffusion separately is that, as will be seen, the inadequacies of a proposed computational method can be more easily detected and eliminated. If the procedure is to be useful, a large number (k) of successive alternate rounds of simulated sedimentation and diffusion should reproduce the actual physical situation at a time k At after the establishment of the specified initial conditions. The effect of the use of a discontinuous model may be tested in two ways. First, a relatively simple case may be selected for which an analytic solution is available. A concentration distribution or boundary shape is calculated using the discontinuous model to be tested. The result is then compared with the exact solution. For the more
249
250
DAVID
interesting cases where the correct result is not known, several separate model calculations can be done, dividing the total time of the experiment to be simulated into different numbers of time intervals; the impact of the discontinuous nature of the description can be evaluated by determining to what degree the results are altered by the arbitrary selection of the number of steps into which the total process is divided. During an ultracentrifuge run, a sedimenting boundary spreads as a consequence of the diffusion of the sedimenting component. In order to simulate such an experiment, it is therefore necessary to develop a computational model for diffusion in the sector-shaped cells generally used with the ultracentrifuge. Vink (3) has reported a simulation procedure for concentrationindependent free diffusion in one dimension in a rectangular cell. In the present paper, the accuracy of Vink’s model is evaluated by comparing the results obtained by its use with the exact solution of Fick’s second law. The method is then suitably modified to deal with cases where diffusion occurs in a sector-shaped cell or where the diffusion coefficient varies with concentration.
J. COX
continuous
analog of Fick’s first law:
SIMULATION AND
PROCEDURES RESULTS
Concentration-independent diffusion in a rectangular cell. The simulation procedure for this simplest case has been published by Vink (3). The description of the method is repeated in detail here to make intelligible the modified procedures described later in the paper. The rectangular cell, with crosssection A, is divided into n compartments by n + 1 boundaries equally spaced along the long axis of the cell; thus, box i is enclosed by boundaries i and i + 1. The interval between the boundaries is Ax, and the distance between the midpoints of adjacent boxes is also Ax. The flow of mass into box i from box i - 1 (AmI) during a short time interval At is given by a dis-
Ci) )
(la)
where D is the diffusion coefficient and Ci is the concentration in compartment i. The mass flow into box i from box i + 1 (Amz) is given by the similar expression:
Am;= D.A.&
(“+’ - “) . (lb) Ax
The sum of Am, and Am, is the change, during At, in the mass of the diffusing component in box i. Normally, one of the Am terms is positive and the other is negative; the sum of Am, and Am2 may be either posi-tive or negative. The change in the concentration of the diffusing component in box i is: ACi
=
AmI + Amz = D.At. A.Ax
;
(&+I - Ci) _ (Ci - ci-1) Ax
-
Ax:
1
@a)
Equation (2a) will be recognized as a discontinuous “translation” of Fick’s second law :
dc -=Dol”c
METHODS Programs incorporating the simulation procedures to be described were written in FORTRAN63 and the calculations were done by a Control Data Corporation 1604 digital computer.
a-‘,,
= D.A.At
Anal
dt
Equation putational ACi
(2b)
dx2 *
(2a) can be simplified purposes: =
,-~(Ci+l
+
ci-1
-
for com-
2ci)
D-At ff=(az)z* Since there can be no mass flow across boundaries 1 and n + 1 at the ends of the cell, the appropriate expressions for compartments 1 and n are: AC1 = a(& - C,)
(24
AC,
(24
= LY(&~
-
c,).
A diffusion experiment is simulated by the following procedure. Given the diffusion coefficient, a suitable box dimension and time interval are selected, and a! is calculated. The initial concentration distribution is specified; for example, if there is initially a
SIMULATION
sharp boundary between solvent and solution at the center of a lOO-box system, the distribution is described by: c; = 0, i = 1;**50 (3)
ci = co, i = 51;.*100,
where Co is the initial concentrat.ion on the solution side of the boundary. The procedure is not restricted to t.his one initial boundary shape. With the initial concentrat.ions specified, Eqs. (2c, 2d, and 2e) arc applied to each compartment. When the calculation has been done for each box, a new set of conccntrations ir calculated. Ci’ = Ci + ACi.
(4)
The set of C’i values defines the concentration distribut.ion after time At. The new concent,rat,ions arc insert.ed into Eqs. (2c, 2d, and 2e), and the calculation is repeated to give the distribution after time 2eAt. Afler h: such transfers, the calculated concentration distribution should reproduce that which would have been observed in a real experiment. k At aft.er the establishment of the initial conditions. The results may bc conveniently converted to a concent,ration gradient, distribution :
0
dc (Ci - Ci-1) d; i = Ax
’
251
OF DIFFUSION
(5)
The gradient subscript in Eq. 5 refers to boundary i; the concentration subscript indicates box i. Gradients are defined only for boundaries 2 through n. Vink (3) has reported using the above procedure for several t.ransfers with a single value of cr; the results appeared promising, in that the calculated boundaries were symand at. Icast approximately metrical Gaussian. The practical use and further elaborat,ion of the method require, however, a detailed examination of the accuracy of the results. In particular, it is essential to know whether the choice of (Y is completely unrestricted. It must be known whether, with a given D, the box interval Ax and the time per t.ransfer At can bc assigned any arbitrary values. It can be shown immediately that. cy is
restricted to values below 0.5. This restricbion is readily understood by considering what happens during the ‘first. t.ransfer in a system with a sharp initial boundary between solvent in boxes 1 through i and solution in boxes i + 1 through n. Examination of Eq. (2~) reveals t.hat, after the first transfer with (Y greater than 0.5, the concentration in box i will be greater t.han that in box i + 1. That is, the concent.ration gradient at the original boundary will be reversed in sign, and a peak in the conccntrat,ion distribution will result. If CYis greater than 1.0, then, after the first transfer, the calculated concentration in box i + 1 will be negative. Thus, with a given diffusion coefficient, physically real results certainly cannot be obtained unless Ax and At arc so chosen that DAM/* is less than 0.5. It remains to see how much effect variation of (Y within the range below 0.5 has on the calculated distribution. A given system may be divided into a large number of narrow compart,ments or a small number of wide boxes; or a given total time Ii A1 may be described by a large k and a small At or vice versa. The fundamental question is whct.her the physical process is accurately depicted with any values of a! and hi. This question was approached empirically by applying the procedure to the diffusion of a component wit.h D = 8.0 X lo-’ cm* per second in a system of 100 boxes, each 0.01 cm wide. At zero time there was a sharp boundary between the solvent and a solution whose concentration was 1.0 mg per milliliter. The gradient curve after 5000 seconds was calculated for this system. Since &
= k $f2
= ka = 40.0,
(6)
the simulation calculation was done by using four different pairs of k and (Y, each selected so that the product of the two was 40.0. In addition, t.hc correct shape of the gradient was calculated by using the exact solut.ion of the appropriate continuity cquation : dc -= dx
co exp (-x2/4Dt). 2(7rDL)“2
(7)
252
DAVID TABLE
SIMULATED
DIFFUSION
I
IN A RECTANGULAR
Diffusion parameter:
0.025
0.05
0.10
0.20
Nzber transfers:
1600
800
400
200
11%- %I1 (cm) 0.00
0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20
CELL
Exact Solution
Concentration gradient (mg/cm~)
4.47 4.36 4.04 3.56 2.99 2.38 1.81 1.31 0.90 0.59 0.37
4.47 4.35 4.04 3.56 2.99 2.38 1.81 1.31 0.90 0.59 0.37
4.46 4.35 4.04 3.56 2.99 2.39 1.81 1.31 0.90 0.59 0.37
4.46 4.35 4.04 3.56 2.99 2.39 1.81 1.31 0.90 0.59 0.37
4.46 4.35 4.04 3.56 2.99 2.39 1.81 1.31 0.90 0.59 0.37
The calculated gradient distributions are displayed in Table I; the results are clear and somewhat surprising. Variation of (Y over a eightfold range from 0.025 to 0.20 has essentially no effect on the shape of the gradient curves. Moreover, the results of the simulation technique are in excellent agreement with the exact solution for all values of ar; the small differences observed are within the range of errors normally encountered in reading t,he schlieren record of a diffusion experiment. The conclusion may be drawn that any convenient choice of Ax and At will be acceptable for the simulation of diffusion by the Vink procedure if the resulting value of (Y is between 0.025 and 0.2; the choice may in fact be more flexible than this. Variation of (Y within me calculation. It is now necessary to modify the simulation procedure described above in order to make it suitable for describing the two cases of immediate interest : concentration-dependent diffusion and diffusion in a sectorshaped cell. Although the specific procedures to be used in these two cases are quite different, the two have one feature in common: the diffusion parameter a! must vary from point to point in the cell. In the sector cell, the cross sections of the successive boundaries and the volumes of the successive boxes increase outward from the axis of the sector. When the diffusion coefficient depends on the local concentration of the solu-
J. COX
tion, and the concentration changes across the cell, then D and thus (Y must vary from point to point. Before proceeding to deal with either case, it is necessary to find out whether the use of different diffusion parameters in different parts of the cell will, of itself, deform the calculated boundary. The simplest way to test this possibility is to recast the rectangular-cell, constant-D case in such a way that (Yvaries from one compartment to the next; this can be done by dividing the cell into compartments of different instead of equal widths. A series of n + 1 boundary positions is specified; an expression analogous to Eq. (2a) is derived for the change, during the time interval At, of the concentration in box i, enclosed between boundaries i and i + 1: DAt Aci = (Xi+1 - Xi)
(8) (Zi - f&l) (is+1 - 5%) 1.
cc,+1 - CJ [
(Ci - ci-1)
In Eq. (8), xi is the position of boundary i; Zi is the position of the midpoint of box i. It should be stressed that the assignment of different widths to the compartments is entirely artificial; the physical process being simulated is still concentration-independent diffusion in a rectangular cell. The curve plotted from the results should therefore be identical to one obtained using a constant box interval; in particular, the calculated gradient curve should be symmetrical and Gaussian. Table II shows the result of one such test. A lOO-box cell was constructed as follows: the width of box 1 was taken as 0.01 cm. The widths of the other boxes were specified by:
(xi+1 - Xi) = (Xi - Xi-l) . g.
(9)
Thus, the width of box i was 0.01 (6.01/ 6.00)i-1, and the boxes increased progressively in width from 0.0100 to 0.0118 cm. Since (Y is inversely proportional to (Ax>~, a! at box 1 was about 1.4 times its value at box 100. The simulation procedure of Eq. 8 was applied for 1000 transfers with DAt equal to 5 X lo+. The resulting gradient distribution was tested for sym-
SIMULATION TABLE
II
TABLE
DIFFUSION SIMULATION WITH CONTINUOUS VARIATION OF Box DIMENSION, TEST FOR SYMMETRY OF GRADIENT CURVE” Box boundaries in lower limb lb - all (4
Concentration gradients in lower limb (mg/cm4) Observed
0.00000 0.01085 0.02170
3.994 3.970 3.900
0.03260 0.04355 0.05450 0.06545 0.07640 0.08735 0.09835
3.786 3.631 3.440 3.221 2.979 2.722 2.458
0.10940 0.12045 0.13150 0.14260 0.15370 0.16480 0.17595 0.18710 0.19825 0.20945 0.22065
2.191 1.929 1.679 1.442 1.224 1.026 0.849 0.694 0.560 0.446 0.350 0.272 0.208 0.158 0.117 0.087 0.063 0.045 0.032 0.023 0.016
OF DIFFUSION
Calculated by interpolation in upper limb
3.994 3.970 3.900 3.786 3.631 3.441 3.221 2.980 2.722 2.457 2.190 1.929 1.677 1.441 1.222 1.024 0.847 0.692 0.558 0.444 0.349 0.271 0.208 0.157 0.117 0.087 0.063 0.045 0.032 0.023 0.016
DISCONTINUOUS Box number 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
III
CEANQE IN Box
DIMENSION
Concentration @g/ml) Variable AZ
Constant Ax
1.000 0.999 0.999 0.998 0.997 0.996 0.993 0.989 0.984 0.977 0.966 0.953 0.935
1.000 0.999 0.999 0.998 0.997 0.995 0.993 0.989 0.984 0.976 0.966 0.952 0.934
0.912
0.911
0.8&l 0.850 0.809 0.763 0.711 0.655 0.594 0.532
0.883 0.849 0.808 0.762 0.711 0.654 0.594 0.532
this case, the box dimension was 0.01 cm for all compartments on the solution side of the 0.24315 initial boundary and 0.02 cm for all compart0.25440 ments on the solvent side. There was a 0.26570 sudden doubling of the box dimension and, 0.27700 thus, a fourfold change in (Y at the center of 0.28830 the diffusing boundary. The concentration 0.29965 0.31100 distribution on the solution side of the 0.32235 boundary after 200 transfers with D At 0.33375 = 1.0 X 10m6is shown in Table III; the initial concentration was 1.0 mg per millidifa Simulated concentration-independent liter. The distribution calculated with the fusion in a rectangular cell. Compartment disame parameters in a set of compartments mensions selected as described in the text. 1000 0.01 cm wide on both sides of the boundary transfers with DAt = 5 X 10-B. is shown for comparison. It is apparent t’hat the correct result is obtained for the recmetry by reflecting the upper limb of the curve through the original boundary posi- tangular-cell, constant-D case, even when (Y is arbitrarily assigned different values in diftion and examining the result to see if it ferent parts of the cell. Simulation procecoincided with the lower limb of the observed curve. As Table II indicates, the dures may now be developed for diffusion in gradient curve was almost perfectly sym- a sector cell and for concentration-dependent diffusion. The results of Tables II and metrical. The sensitivity of the rectangularcell, constant-D simulation to variation in cr III lend some confidence that the distribuwas subjected to a more rigorous test by de- tions calculated for these cases will not be distorted by the fact that the diffusion pafining still another set of compartments. In 0.23190
254
DAVID
rameter happens to vary with position in the cell. Diflusion in a sector-shaped cell. The physical situation to be described is the spreading of a boundary in a synthetic boundary cell in an ultracentrifuge operated at low speed. In the usual notation, the sector angle is 8, and the height of the sector is b. The cell is divided into n compartments by n + 1 cylindrical surfaces located at various distances r from the axis of the sector; boundary 1 is the meniscus, and boundary n + 1 is the cell bottom. The radial dimensions of the various compartments need not be the same. The area of boundary i is :
(10)
Ai = mi.
The volume of box i, between boundaries i and i + 1, is: vi = W&i+1
- ri),
(11)
where pi is the radial position of the midpoint of box i: i;,
=
b-i
z
+
ri+11
2
(12)
*
During the time At, the flow of mass into box i from box i - 1 across boundary i is:
and the flow of the diffusing component into box i from box i + 1 across boundary i+ 1 is: (Cc+1 - Cd Am2 = D. At. bf?ri+l (pi+l _ pi) .
(13b)
The change in the concentration of the diffusing component in box i during time interval At is: AC,
z-
_
Am1 + Am2 Vi
=
J. COX TABLE IV COMPARISON OF SIMULATED DIFFUSION RECTANGULAR AND SECTOR-SHAPED CELLS
Concentration gradient (mgjcm’)
Since the box dimensions during the calculation, the putations may be reduced diffusion parameters bi and
1
(14)
do not change number of comby defining the fi:
Rectangular cell
sector cell
0.508 0.625 0.762 0.918 1.094 1.290 1.504 1.735 1.980 2.234 2.493 2.750 3.000 3.235 3.449 3.636 3.788 3.901 3.970 3.994 3.970 3.900 3.786 3.631 3.441 3.221 2.980 2.722 2.457 2.190 1.929 1.677 1.441 1.222 1.023 0.847 0.692 0.558 0.444
0.516 0.634 0.773 0.930 1.107 1.305 1.521 1.753 1.998 2.253 2.511 2.768 3.017 3.251 3.464 3.647 3.797 3.907 3.973 3.993 3.967 3.893 3.776 3.619 3.427 3.206 2.962 2.704 2.439 2.172 1.911 1.660 1.425 1.207 1.011 0.835 0.682 0.549 0.438
6.3126 6.3231 6.3337 6.3442 6.3548 6.3654 6.3760 6.3866 6.3973 6.4079 6.4186 6.4293 6.4400 6.4508 6.4615 6.4723 6.4831 6.4939 6.5047 6.5155 6.5264 6.5373 6.5482 6.5591 6.5700 6.5810 6.5919 6.6029 6.6139 6.6250 6.6360 6.6471 6.6581 6.6692 6.6804 6.6915 6.7026 6.7138 6.7250
DAt c;i(ri+1 - f-i>
_ ri+1cci+1- Ci> (Pg.1 - Pi)
IN
bi = I
Ji =
DAt ri+l Pi(Ti+1 - ri)(fi+l _ Ti(Ti+l
thus simplifying ACi = fi(C,,
- Pi) (15)
DAt r; -
ri)(Pi
-
f&l)
’
Eq. (14) to give: - Ct> + b;(Ci+l -
Ci). (16)
SIMULATION
The condition that no mass flows across boundaries 1 and n + 1 at the ends of the cell is: fl = b,+1 = 0.
255
OF DIFFUSION
(17)
Diffusion in a sector-shaped cell is simulated in the following way. A convenient set of boundary locations ri is specified. The locations of the midpoints of the boxes, pi, are calculated. A time interva.1 At is selected and, with D given, the diffusion parameters bi and fi are determined for each box; At must be so chosen that each of the diffusion parameters is less than 0.5, and preferably lies within the range 0.025-0.2. The initial concentration distribution is specified, and Eq. (16) is applied to each box. The concentration distribution at time At is calculated, and the process is repeated k times. The calculated concentration distribution should reproduce the physical situation k At after the beginning of a real experiment. Table IV shows gradient distribution resulting from such a calculation. The set of boundaries chosen for the simulation happened to be the ones specified by Eq. (9) ; i.e., the width of the boxes varied across the cell. Thus the results can be compared directly with those shown in Table II. There is, of course, no necessity for specifying a variable cell dimension in order to describe the sector-cell case. The model cell was placed between 6.00 and 7.08 cm from the axis of the sector; the model was thus very much like the usual ultracentrifuge cell. As in the experiment of Table II, the initial sharp boundary was placed between the 50th and 51st boxes; the initial concentration on the solution side of the boundry was; 1.0 mg per milliliter, and diffusion was simulated by 1000 transfers with D At equal to 5 X 10-6. Inspection of Table IV reveals that the effect ot the sector shape of the cell is small but noticeable. The shape ot the gradient curve is not detectably changed, but the whole curve is shifted very slightly (about 7~) toward the narrow end of the sector. Concentr&on-dependent diffusion in a rectangular cell. In simulating this important case, the diffusion coe@icient at infinite dilution is specsed, and the dependence of D
on concentration series :
is described by a power
D = DO (1 + k# + kzC2 +. . . -).
(18)
The concentration which determines the local diffusion coefficient at boundary i is taken to be the mean of the concentrations on the two sides of the boundary, in boxes i and i - 1. Diffusion where D is a linear function of concentration has been considered in detail; in this case, the local diffusion coefficient at boundary i is: Di = Do 1 + k (ci +2 “-‘))
.
(19)
It is convenient to divide the cell into n compartments by n + 1 equally spaced boundaries, as in the concentration-independent case. Equation (2a) is modified to yield a working equation suitable for concentration-dependent diffusion: 1+
ACi = 00
k(Ci
+
Ci-1)
>
K
* (Ci-1 - C;) +
1 + ,“(” * (&+I
: “‘+l)> -
Ci>
(2Oa)
1 7
where
(20b) The fact that no mass crosses the ends of the cell, boundaries 1 and n + 1, is stated as follows : k(Cl + 2
c2,>
(2Oc)
* cc2 - Cd G’Od)
*(C?r,- cn>. A concentration-dependent diffusion experiment was simulated with the working Eqs. (20a), (~OC), and (20d); the result is shown in Fig. 1. The model cell consisted of 100 boxes 0.01 cm wide. The sharp boundary
256
DAVID
J. COX
50.0 7 E 9 40.0 If ‘: .? 30.0 ‘0 E!
W
.z 20.0 ‘0 L ‘; 0” 6 10.0 0
-3.0
-2.5
-2.0
-1.5
-1.0
-0.5
0
X-X,
to.5
t1.0
t1.5
t2.0
+2.5
+3.0
3.5
(mm)
FIG. 1. Simulation of diffusion with D a linear function of concentration. Diffusion coefficient varies from 2.5 X lO+ cm2/second at C = 0 to 7.5 X 10-7 cm2/second at C = 10 mg/ml. Triangles show the concentration gradient distribution after 8000 seconds. Circles show the gradient distribution resulting from the diffusion of a component with D = 5.0 X 10-7 cmz/second for 8000 seconds.
between the solvent and a solution at 10 mg per milliliter was formed between boxes 50 and 51. For illustrative purposes, an extreme case was selected in which the diffusion coefficient varied linearly with concentration from 2.5 X lo-’ cm2 per second at infinite dilution to 7.5 X lO-’ cm2 per second at 10 mg per milliliter. The time interval per transfer was 10 seconds, and the simulation was carried out for 800 transfers. The diffusion for 8000 seconds of a component with D equal to 5.0 X lo-’ cm2 per second at all concentrations was also simulated; the resulting gradient curve is shown in Fig. 1 for comparison with the analogous concentration-dependent experiment. Gosting and Fujita (4) have reported an exact solut’ion of the continuity equation for the case where the diffusion coefficient is a function of concentration. The expression obtained for the concentration gradient distribution contains a number of terms in the form of a power series in initial solute concentration. Costing (5) has shown a plot based on this solution which is analogous to Fig. 1; it compares the gradient curves for concentration-dependent and concentrationindependent diffusion. However, the plot for the concentration-dependent case is an approximation; it takes account of only the
TABLE COMPUTER
V
TIME REQUIRED FOR VARIOUS SIMULATION PROCEDURES
Process simulated
NO. transfers
Rectangular-cell, constant D Sector-cell, constant D Rectangular cell D - Do(1 + kc) Rectangular cell D - Do(1 + MY)
1000
1 minute
14 seconds
1000
1 minute
20 seconds
1000
1 minute
35 seconds
500
3 minutes
17 seconds
Computer time
first two terms of the exact solution. Fig. 1 and Gosting’s Fig. 6 are quite similar. The gradient distribution for the concentrationdependent system is markedly skewed toward the region of the lower diffusion coefficients. The plot for the concentrationdependent case crosses the analogous curve for the system with a constant diffusion coefficient at three points. In the graph presented by Gosting, the crossover points occur at the peak of the symmetrical concentration-independent plot and at two other points equidistant from the maximum of the constant-D gradient distribution. In Fig. 1, the middle crossover point is dis-
SIMULATION
placed from the peak of the Gaussian constant-D plot toward the region of lower concentrations, and the other two crossover points are not symmetrically arranged around the peak. The differences are probably due to the neglect of the higher terms of the analytic solution in the construction of Gosting’s Fig. 6. The sum of these higher terms is apparently negative and approxirnately symmetrical (4). No great amount of computer time is needed for calculations of this kind. Table V shows the time required to process simulation programs of several types. DISCUSSION
The procedures described in this paper were devised as a part of a more general project leading to the development of methods for the simulation of sedimentation in the ultracentrifuge. In this context, the major results of the work are that the Vink procedure describes diffusion very accurately, and that the method can easily be modified to deal with diffusion in a typical ultracentrifuge cell. The simulation of diffusion may be of some practical use in itself. Exact solutions of Fick’s second law are, of course, available and are widely used in systems where the variation of the diffusion coefficient with concentration is negligible; the simulation procedure has no important practical use in these simple cases. Where the diffusion coefficient varies with concentration, however, the methods presently in use leave something to be desired. Exact solutions are available for handling the results of diffusion experiments when D is a function of solute concentration (4). These exact solutions are somewhat cumbersome in practice, however; their use involves a great deal of laborious computation. As a consequence, methods for handling such data generally rely on approximations to the exact solution. The most commonly used simplification consists of discarding the higher terms in the series solution (6). The total effect of such a mathematical approximation on the accuracy of the calculated results is generally not obvious, but one assumption implicit in this procedure is that the dependence of D
OF DIFFUSION
257
on powers of concentration above the first is negligible. When diffusion coefficients at different solute concentrations are calculated using the standard approximate treatments, it is not unusual to find that D does not vary linearly with concentration. Infinite-dilution diffusion coefficients derived by extrapolation from such data are of uncertain accuracy, particularly if the concentration dependence of the diffusion coefficients happens to be steep. The simulation technique may be suggested as a method which would place the calculation of diffusion coefficients for concentration-dependent systems on firmer ground. An infinite-dilution diffusion coefficient and the parameters describing the dependence of D on concentration might be calculated by the usual approximate methods. A diffusion experiment could then be simulated by using these approximate parameters and the result of the simulation could be compared with the observed boundary shape. If the calculated and observed boundaries did not coincide, the diffusion parameters could be adjusted to give a better fit. In principle, there is no reason why the computer could not be presented with the experimental concent,ration or gradient distribution and instructed to find a Do and regression coefficients that would generate the observed results. However, a program of this type has not yet been written. The most attractive feature of the simulation technique is its great flexibility. The exact solution of the continuity equation is very much complicated by a two- or threeterm concentration dependence of the diffusion coefficient, but little more computer time is required for the description of such a system by the simulation technique t,han for the simple constant-D calculation. Moreover, the initial boundary in a diffusion experiment is usually not perfectly sharp. The methods generally used for diffusion calculations formally handle a blurred initial boundary in terms of a time error. This procedure contains an implicit assumption about the shape of the initial boundary; for example, it assumes, for the constant-D case, that the initial boundary is Gaussian.
258
DAVID J. COX
The application of t.he simulation technique requires only that. the shape of t.he initial boundary be known, so that it may be specified at the beginning of the calculation. REFERENCES Theory of Sedi1. FUJITA, II., “Mathematical mentation Analysis.” New York, Academic Press (1962).
2.
J. L., AND KEGELES, G., J. Phjys. Chem. 66, 1761 (1961). 3. \'Ih'K, II., nCk& Chem. &and. 18, 409 (1964). 4. GOSTIN, L. J., ASD FLTJIT.~,I-I., J. Am. Chem. sot. 79.1359 (1957). 5. (doSTISG, 1~. J., A&an. Protein Chem. 11, 429 (1956). 6. CHEETH, J. M., J. Am. Chem. Sot. 77. 64% (1955). BE?.IIUNE,