A Precise Boundary Element Method for Macromolecular Transport Properties SERGIO ARAGON
Department of Chemistry & Biochemistry, San Francisco State University, 1600 Holloway Ave., San Francisco, California 94132 Received 7 October 2003; Accepted 9 February 2004
Abstract: A very precise boundary element numerical solution of the exact formulation of the hydrodynamic resistance problem with stick boundary conditions is presented. BEST, the Fortran 77 program developed for this purpose, computes the full transport tensors in the center of resistance or the center of diffusion for an arbitrarily shaped rigid body, including rotation-translation coupling. The input for this program is a triangulation of the solvent-defined surface of the molecule of interest, given by Connolly’s MSROLL or other suitable triangulator. The triangulation is prepared for BEST by COALESCE, a program that allows user control over the quality and number of triangles to describe the surface. High numerical precision is assured by effectively exact integration of the Oseen tensor over triangular surface elements, and by scaling the hydrodynamic computation to the precise surface area of the molecule. Efficiency of computation is achieved by the use of public domain LAPACK routines that call BLAS Level 3 hardware-optimized subroutines available for most processors. A protein computation can be done in less than 10 min of CPU time in a modern Pentium IV processor. The present work includes a complete analysis of the sources of error in the numerical work and techniques to eliminate these errors. The operation of BEST is illustrated with applications to ellipsoids of revolution, and Lysozyme, a small protein. The typical numerical accuracy achieved is 0.05% compared to analytical theory. The numerical precision for a protein is better than 1%, much better than experimental errors in these quantities, and more than 10 times better than traditional bead-based methods. © 2004 Wiley Periodicals, Inc.
J Comput Chem 25: 1191–1205, 2004
Key words: boundary element; hydrodynamic resistance problem; rotation-translation coupling
Introduction The structure of molecules in the solid state can be readily determined by X-ray diffraction methods,1 even for very large macromolecules, as long as they can be crystallized. In most instances, however, the molecules are of interest because of their function while immersed in liquid aqueous environment. The size, shape, and available conformations of molecules, especially those that can form supramolecular aggregates, polypeptides, or multisubunit proteins, or nucleic acids, may be different in solution compared to a crystalline state.2 Thus, many methods exist to probe the structure and dynamics of molecules in solution, including dynamic laser light scattering,3 transient electric birefringence,4 fluorescence polarization anisotropy,5 fluorescence photobleaching recovery,6 electron spin resonance,7 and nuclear magnetic resonance.8 These modern techniques emphasize the dynamics of the molecules in solution at time scales characteristic of each technique. Classical techniques such as centrifugation and electrophoresis are also important for biomolecules.9 The simplest (and slowest) of the dynamical motions available to molecules in solu-
tion involve translation and rotation. It is these motions that correlate most directly with the overall molecular size and shape. Flexing, conformational changes, and local segmental motions are typically much faster processes,10,11 and they contribute to the shape in an average way. The characteristic time constants for the slower process are often quite well described as diffusive in origin, and thus directly relate to the hydrodynamic friction tensors of the molecule in question. Thus, a hydrodynamic description of macromolecules is of interest in the interpretation of a wide class of experimental measurements. For macromolecules, consideration of the solvent as a continuum is an excellent approximation (with modification of the boundary conditions, it is also excellent for small molecules). Thus, the governing equations for the computation of the hydrodynamic transport properties are the Navier–Stokes equations of fluid flow. In the limit of small Reynolds number, as appropriate Correspondence to: S. Aragon; e-mail:
[email protected] Contract/grant sponsor: National Institutes of Health (MBRS SCORE Program); contract/grant number: S06 GM52588.
© 2004 Wiley Periodicals, Inc.
1192
Aragon
•
Vol. 25, No. 9
•
Journal of Computational Chemistry
for the diffusion process, the equations are known as the Stokes or creeping flow equations.12 The spherical case, along with spheroidal shapes, rings, and polygons, are among the few shapes for which the creeping flow equations of hydrodynamics can be solved exactly.13 Most molecules do not have such high symmetry or smooth contours so that a method that takes the precise molecular shape into account is needed. Methods based on a set of hydrodynamically interacting beads were pioneered by Bloomfield14 and coworkers in the late 1960s. The work based on bead models progressed rapidly in the following decade, and it has been reviewed by Garcia de la Torre and Bloomfield.15 The idea is simple and appealing. Due to the intrinsic limitations of the method, most workers followed Bloomfield’s example and represented whole functional groups in nucleic acids or amino acid residues in proteins as a single frictional bead.16 The beads were located along the backbone of the polymers, and thus the general contour of the biomolecule was generated in a coarse manner. Even though this method sounds crude, it achieved a reasonable measure of success in extracting shape and size information from the more classical intrinsic viscosity and sedimentation coefficient measurements.17 In bead models, the hydrodynamic interaction (HI) between the beads is the crucial ingredient in the computation of the measurable diffusion coefficients. Several authors have obtained good representations of the HI, including dependence on actual bead sizes.15,18 –20 For the crucial case of equal size overlapping beads, the HI tensor was obtained by Rotne and Prager,20 but none is as yet available for the case of overlapping beads of unequal size. In a coarse-grained model, one can either avoid bead overlaps, or select overlapping beads of equal diameter when necessary. When the modeler has at his/her disposal the size of the beads to use in a coarse-grained model, then it is always possible to fit a given experimental measurement. However, it is generally the case that if the model parameters have been set by fitting to a translational diffusion coefficient, then the rotational diffusion coefficient is not accurately predicted, or vice versa.16,21–23 This is an intrinsic limitation of coarse bead models. Ideally, one would like to be able to predict the hydrodynamic transport coefficients to high accuracy knowing the atomic coordinates of the molecule, whether the molecule is small or large. Assigning a Van der Waals radius to each atom in the molecule requires that the “atomic beads” overlap, thus making the computation of the hydrodynamic interaction in a realistic atomistic model impossible. A nonrigorous way to get around this problem has been proposed.24 Pastor and Karplus25 have proposed an alternative atomistic hydrodynamic modeling method. To avoid the overlap problem, they introduced the accessible surface area (ASA) model where the frictional radius of each atom is calculated by measuring its contribution to the molecular Van der Waals surface that is generated by all the atoms in the molecule. Atoms hidden inside the molecule contribute nothing. This type of calculation still depends on the HI tensor, but avoids the overlap problem because the atomic radius is reduced by the ASA model. This method has been used successfully to model small molecules, proteins, and viruses. However, the method is a hybrid implementation of “stick” and “slip” boundary conditions and is thus not a rigorous hydrodynamic treatment. When large proteins are attacked by this method,26 good agreement can be obtained with experimental transla-
tional diffusion coefficients, but the amount of water of hydration can always be adjusted to reach that agreement. With small molecules,25 however, the agreement is actually not that good. The typical accuracy obtainable by the above bead methods varies from 5 to 20%. Greater accuracy than this is required to make an unequivocal statement about the degree of hydration, the solution structure, or the conformation of molecules in solution. Recently, however, Garcia de la Torre et al.27 have further developed the “shell” model originally suggested by Bloomfield. In this method, very small equal sized beadlets are placed on the surface of the body to be modeled, avoiding any overlaps. Then, by assigning multiple beads to the solvent accessible surface of atoms, a bead model could produce an acceptable atomistic transport coefficient in the limit of zero bead size. Because the HI tensors are only approximate, however, this method requires a different parametrization depending on which transport property is being addressed. If we are to focus on the surface of the body, then we are better off using an exact representation of the transport properties without reference to beads at all. This is the objective of this work. Whereas bead methods aim to solve a mobility problem, which cannot be formulated exactly, an alternative method is to solve a resistance problem, which can be formulated exactly. As is shown below, once one has precise friction tensors, it is straightforward to compute the diffusion tensors. In the mid-70s Youngren and Acrivos28 presented an effective method for the numerical solution to the exact surface integral representation of the velocity field for the creeping flow equations. This expression was used to calculate stick-and-slip hydrodynamic transport coefficients for spheroids to high accuracy. The method lends itself easily to numerical analysis as stresses distributed over the wetted surface of the body. The method developed by Hu and Zwanzig29 is a special case of the general formulation of Youngren and Acrivos (YA). A decade later, Wegener30 made a clear exposition of the advantages of this model and advocated its use. The formulation is general, and thus applies to either stick or slip boundary conditions. It is the only method that rigorously includes the slip boundary condition case. In the intervening time, hardly any work appeared in the literature utilizing the YA method. Recently, Allison has presented calculations on the electrophoresis problem based on a simple implementation of the YA method,31–34 and the transport properties of simple shapes using the “slip” boundary condition.35 The YA method has not been used extensively so far for several reasons. First, the integral equation approach, compared to the simplicity of beads, is more mathematically obscure (despite the fact that the basic equations date from the early 1930s and the full modern exposition of their application to hydrodynamics by Ladyzhenskaya36 predates the YA work by more than a decade). The bead methods introduced by Kirkwood37 in the 1940s to describe polymers have had decades of historic success, and have been well entrenched in the field, well before the advent of computers. Bead methods allow simplified calculations that can be carried out by hand, yet yield very useful results. Integral equation methods, on the other hand, need the powerful computers of today for their successful implementation. Today, supercomputer power is available at the individual desktop—the day for the widespread application of the YA method has arrived. We intend to show that the YA method is the most general hydrodynamic method and one that
Boundary Element Method for Macromolecular Transport Properties
is capable of producing unprecedented accuracy in the computation of transport properties of macromolecules that surpasses the accuracy obtained from the bead methods.
In the following we review the equations of the Youngren–Acrivos28 method. For the case of macromolecules, “stick” boundary conditions are appropriate. In this case, the velocity field of the flow, v(y) at position y in the fluid, can be written as an integral over the particle surface (SP),
冕
冘 N
v共yk 兲 ⫽
7
Gkj 䡠 fj
(4)
j⫽1
The centerpiece of this set of equations is a set of N completely known 3 ⫻ 3 matrices of coefficients that contain all geometric information, the integrals of the Oseen tensor over a surface patch,
Theory: The Boundary Element Method for Stick Boundary Conditions
v共y兲 ⫽ uo 共y兲 ⫹
1193
7
T共x, y兲 䡠 f共x兲dSx
(1)
sp
where uo (y) is the flow velocity of the fluid if the particle was not there (which can be taken to be zero for diffusive motion), and T(x, y) is the Oseen hydrodynamic interaction tensor. The surface stress force, f(x) is the unknown quantity that we must obtain. Once this quantity is known, the transport properties of the macromolecule can be directly computed, as shown below. The Oseen tensor, is given by,12,13
7
G kj ⫽
冕
7
T共x, yk 兲dSx
(5)
⌬j
The set of 3N equations can be written all at once,
冤冥 冤 v1 ··· ··· vN
7
G11 ··· ⫽ ···
3Nx1
7
GN1
7 ··· ··· G 1N ··· ··· ··· ··· ··· ··· 7 ··· ··· G NN
冥 冤冥 f1 ··· ··· fN
3Nx3N
(6) 3Nx1
7
7
T共x, y兲 ⫽
冋
册
1 7 共x ⫺ y兲共x ⫺ y兲 I⫹ 8兩x ⫺ y兩 兩x ⫺ y兩2
(2)
It is important to note that in bead hydrodynamics, the Oseen tensor is only the first term in an infinite series expansion of the interaction between two beads centered at x and y, respectively. However, when the hydrodynamics is expressed as a continuous integral over the surface of the body, the tensor is an exact representation of the hydrodynamic interaction of the infinitesimal surface elements. Thus, the starting expressions for the calculation, unlike the bead modeling case, are exact;28,30 moreover, the equation is applicable to an arbitrarily shaped body. Because eq. (1) is an integral equation, the solution requires an approximate numerical method. The method, however, can be iterated to obtain arbitrary precision. The first step is to discretize the surface by replacing it with a collection of N patches that smoothly tile the molecular surface (methods for doing this will be discussed below). Then we can write,
冘 N
SP ⫽
⌬j
(3)
j⫽1
We place the coordinate xj at the center of the small patch ⌬ j and take the surface stress force f(x) to be a constant over the entire patch area. This is the basic approximation: it is clear that it will become a better and better approximation as the patch is made small. Thus, an extrapolation to zero size patch leads to a very precise value for the transport properties. With this approximation, eq. (1) becomes a set of 3N equations for 3N unknowns f(x),
from which the unknown surface stress forces can be readily 7 obtained by matrix inversion of the 3N ⫻ 3N super matrix G, 7
⫺1 关f兴3Nx1 ⫽ 关G兴3Nx3N 关v兴3Nx1
(7)
The total force and torque on the body can be computed from the surface stress forces, and these are directly related to the 7 friction tensors (K) of the body,
冘 N
F⫽
7
7
fj 共x兲⌬j ⫽ ⫺Ktt 䡠 vp ⫺ Ktr 䡠 p
(8)
j⫽1
冘 N
T⫽
7
7
xp ⫻ fj 共x兲⌬j ⫽ ⫺Krt 䡠 vp ⫺ Krr 䡠 p
(9)
j⫽1
The particle can be assumed to have specific translation velocity vp and angular velocity p (for example p ⫽ 0 and vp ⫽ (v x , 0, 0)) to solve the above equations. Thus, six calculations suffice to determine all components of the friction tensors. The friction tensors form part of a larger 6 ⫻ 6 tensor that contains information about the pure translational friction (tt), the pure rotational friction (rr), and the coupling that may exist between these (rt and tr). There are actually only three independent friction tensors because 7 7 the Ktr tensor is the transpose of the Krt tensor. This coupling is insignificant unless the body has a screw-like axis of symmetry.38 The diffusion tensors are finally obtained from the friction tensors by an easy 3 ⫻ 3 matrix inversion, 7
7
7
7
7
7
7
7
7
7
⫺1 D tt ⫽ kT关Ktt ⫺ Ktr 䡠 Krr 䡠 Krt 兴⫺1
⫺1 D rr ⫽ kT关Krr ⫺ Ktr 䡠 K⫺1 tt 䡠 Krt 兴
(10) (11)
1194
Aragon
•
Vol. 25, No. 9
•
Journal of Computational Chemistry
These formulas need one last correction. It so happens that the value of the rotational diffusion tensor in eq. (11) is independent of the choice of origin for the location of the surface patches of the body; however, the same is not true for the translational diffusion tensor in eq. (10). There is one preferred center, called the “Center of Diffusion,” which must be selected to get the correct measurable value for translation. The procedure to achieve this entails an analysis of the transformation properties of the tensors with respect to translations that has been given by Brenner.38 We will omit most of the tensor equations for brevity. Suffice it to say that in the 7 center of diffusion, the Dtr tensor, which is given by, 7
7
7
7
D tr ⫽ ⫺K⫺1 tt 䡠 Ktr 䡠 Drr
(12)
must be symmetric. If the origin that we used is not at the center of diffusion, the tensor in eq. (12) will be asymmetric. The asymmetry present in the rotation-translation coupling tensor of eq. (12) can be used to calculate the translation vector to the center of 7 diffusion to correct the value of Dtt . This procedure has been implemented in our Fortran programs described below.
Implementation of the BE Method: BEST and COALESCE The implementation of a precise computational method must consider two closely tied aspects: the generation of the discretized surface, and the hydrodynamic computation given those surface
Figure 1. Tessellated benzene (1068) facets. [Color figure can be viewed in the online issue, which is available at www.interscience. wiley.com.]
Figure 2. The extrapolation of the perpendicular rotational diffusion coefficient of benzene under stick boundary conditions. The squares correspond to the area corrected computations, while the diamond tesselations have an increasing area as the number of triangles, N, increases.
elements. We begin our discussion with the sources of error in the hydrodynamic computation. There are four fundamental sources of error one must consider in a precise implementation of the boundary element method. These are: (1) the discretization error, (2) the integration error, (3) ill-conditioning and the shape and area error, and (4) the round-off error. The discretization error arises from the approximation of constant surface stress over a surface boundary element. As shown below, this error can be eliminated to achieve accuracies of 0.05% or better by extrapolation to large enough number of surface elements. This error cannot be avoided, but it can be made small, even without extrapolation, if the other errors are minimized or eliminated. Thus, let’s look at these other sources of error in detail. The integration error arises from any approximations adopted in the integration of the Oseen tensor over the triangular elements. To study the integration and discretization errors, it suffices to consider a simple nonspherical shape. Figure 1 shows a typical surface discretization of benzene performed with the programs described below. Benzene has been selected because it is nonspherical, and for purely illustrative purposes. The proper hydrodynamic treatment of benzene requires slip boundary conditions.3 Thus, the numbers we obtain in this “stick” computation will not be compared to experiment. In Figure 2, we have plotted the perpendicular component of the rotational diffusion tensor for “stick” benzene for two instances. The diamonds correspond to a straightforward calculation where the G super matrix of eq. (6) has been computed with insufficient numerical precision. In particular, the diagonal matrix elements have been computed to doubleprecision accuracy, but the off diagonal elements have been computed to only single-precision accuracy. The immediate consequence of this choice is that the plot of the hydrodynamic properties vs. the inverse of the number of triangles is extremely noisy, and the extrapolation to infinite N has unacceptably low precision and accuracy. On the other hand, the data represented by squares in Figure 2, can be extrapolated to very good numerical precision. This data
Boundary Element Method for Macromolecular Transport Properties
1195
Table 1. The Perpendicular and Parallel Tensor Components for Translation and Rotation of Ellipsoids.
p Oblate 2 4 8 30 Prolate 1 1 2 1 4 1 8 1 30
Dtt ⬜
%
Dtt 㛳
%
Drr ⬜
%
Drr 㛳
%
0.84104 0.48853 0.26664 0.076401
0.004 0.004 0.01 0.02
0.73655 0.38441 0.19514 0.052341
0.02 0.02 0.02 0.02
0.22105 0.033933 4.5054 10⫺3 8.7111 10⫺5
0.06 0.07 0.07 0.02
0.17728 0.027774 3.9623 10⫺3 8.3765 10⫺5
0.01 0.03 0.04 0.07
1.33356
0.02
1.33356
0.02
1.0003
0.03
1.0003
0.03
0.96702
0.007
1.10782
0.03
0.3323
0.02
0.61999
0.03
0.64839
0.006
0.83443
0.003
0.073631
0.01
0.34671
0.003
0.40954
0.0007
0.57596
0.04
0.013296
0.008
0.18224
0.04
0.15318
0.01
0.23973
0.06
3.9941 10⫺5
0.02
0.049895
0.13
for prolate ellipsoids, p ⬍ 1
(13e)
The % columns indicate the error by comparison with the analytic formulas.
was computed with the area correction, described in detail below. It is important to note that precision in the computation of the super G matrix is the same in both cases for Figure 2; the “area correction” is what makes the difference. Other workers who have applied the boundary element method to hydrodynamics have either not intended to produce very high precision, or not attempted to extrapolate to very large number of surface elements, so this problem has not been previously explicitly addressed. As Youngren and Acrivos28 admonished, it is worth our while to make the integrations of the G matrix as accurate as possible. Transport Properties of Ellipsoids
To study precision and accuracy further, we turn to the computation of the transport tensors for ellipsoids of revolution and compare to the known analytic formulas. The theoretical expressions for the explicit tensor components are available in Kim and Karilla.13 (Note that the typical Perrin formula for translation gives one third the trace of the tensor, while we compute all the individual components here.) These expressions are functions of the axial ratio p ⫽ b/a, where a is the semiaxis of revolution. We have also omitted a factor of kT/8 for convenience in both analytical and numeric work, and the magnitude of a is in arbitrary units. For prolate ellipsoids, b has been set to 1, and a varied, while for oblate ellipsoids, a has been set to 1 and b varied. 1 G关 p兴共2 ⫺ p 2兲 ⫺ 1 a 1 ⫺ p2
(13a)
1 G关 p兴共2 ⫺ 3p 2兲 ⫹ 1 2a 1 ⫺ p2
(13b)
1 3共1 ⫺ p 2G关 p兴兲 a 3p 2 2共1 ⫺ p 2兲
(13c)
1 3共G关 p兴共2 ⫺ p 2兲 ⫺ 1兲 a3 2共1 ⫺ p 4兲
(13d)
Dtt 㛳 ⫽
Dtt ⬜ ⫽
Drr 㛳 ⫽
Drr ⬜ ⫽
冋
G关 p兴 ⫽ log
1 ⫹ 冑1 ⫺ p2 p
册冒 冑
1 ⫺ p2 ,
G关 p兴 ⫽ ArcTan关 冑p2 ⫺ 1兴/ 冑p2 ⫺ 1, for oblate ellipsoids, p ⬎ 1.
(13f)
Table 1 presents the computation of the parallel and perpendicular components of the translational and rotational diffusion tensors for both prolate and oblate spheroids as a function of axial ratio. The quantities in Table 1 are the result of extrapolations to infinite N on the basis of triangulations performed with Mathematica. Figure 3, shows sample triangulations for an oblate and a prolate spheroid, while Figure 4 shows a typical extrapolation for one of the tensor properties, along with statistical information for the least squares fit. Figure 4 shows that the standard error in the fit is smaller than the actual accuracy; thus, the statistical precision exceeds the accuracy in this calculation. The fit has a small curvature, but it is statistically well defined because the Tstat magnitude is decidedly larger than 1. A linear fit clearly shows (data not shown) nonrandom residuals, and less accuracy in the results. The above fit and its statistics are completely typical of all the computations done on ellipsoids, both prolate and oblate. For each case, four to five computations with varying number of triangles were performed to provide extrapolation data. The only case with a larger error is the needle-like prolate ellipsoid whose axial rotational diffusion coefficient has about twice the error typical of other quantities. This is not unexpected due to the difficulty of triangulating accurately the very fine tip of the needle. The error is barely over 0.1%, even in this worst case. The data in Table 1 and Figure 4 show that, with sufficient precision, the discretization error can indeed be eliminated. The numerical results are essentially exact (with one exception) with errors in the fifth decimal place (0.07% or less), even for high axial
1196
Aragon
•
Vol. 25, No. 9
•
Journal of Computational Chemistry
Figure 3. (a) Triangulations for an oblate ellipsoid of axial ratio 4 with 1104 triangles. Top view on the 1 left and sideways view on the right. (b) Triangulations for a prolate ellipsoid of axial ratio 4 with 1104 triangles. Top view on the right and sideways view on the left. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]
Boundary Element Method for Macromolecular Transport Properties
1197
1
Figure 4. Least-squares fit to Dtt perpendicular vs. 1/N for a prolate ellipsoid with axial ratio 4.
ratios. The computations were done with double precision accuracy for all components of the super G matrix in this case, but still utilizing the “area correction.” The Area Correction: Motivation
It is fairly obvious that the tessellation of the surface yields an area smaller than the true surface area of the ellipsoid, leading to less friction: this is the shape and area error. The data of Table 1 show that this error can also be essentially eliminated with the “area correction,” even with a finite discretization. The area correction is a nonuniform scaling of the super G matrix elements, and is based on the area dependence of the exact integration of the Oseen tensor over a triangle. The details on this scaling are presented below. More importantly, this “area correction” also helps us avoid illconditioned G matrices in this problem. The full study of the issues in ill-conditioned matrices arising in the numerical solution of Fredholm integral equations of the first kind will be addressed in
a separate communication. The area correction also has the advantage that it yields an absolute error that is smaller at every N (the slope of the properties vs. 1/N has decreased), and the graph remains linear. A high precision computation that suffers from the effects of ill-conditioning is illustrated in Table 2 for the case of a prolate ellipsoid of axial ratio 41. The eigenvalues of the translational diffusion tensor are very precise, and the extrapolated values agree with the theoretical values to very high accuracy. However, one of the three components of the rotational diffusion tensor, Drr ⬜ (2), is off by several orders of magnitude while the other two are accurate (the intrinsic viscosity is also unphysical for this case, but this topic will not be discussed here). This phenomenon is observable for both oblate and prolate ellipsoids of axial ratio larger than 1. The sphere does not suffer from this problem (we have not determined how much the axial ratio must differ from unity to observe the ill-conditioning). We have also observed the effects of
1198
Aragon
•
Vol. 25, No. 9
•
Journal of Computational Chemistry
Table 2. An Ill-Conditioned Computation for the Prolate Ellipsoid.
p
Dtt ⬜
Prolate—Ill Conditioned 1 0.64835 4
%
Dtt 㛳
%
Drr ⬜ (1)
%
Drr ⬜ (2)
%
Drr 㛳
%
0.0002
0.8344
0.0002
0.073631
0.009
0.000281
105
0.34658
0.04
The component Drr ⬜ (2) of the rotational diffusion tensor is grossly incorrect. This happens for all the ellipsoid axial ratios studied here, both prolate and oblate. It does not happen for the sphere.
ill-conditioning in computations of other geometrical shapes such as a tetrahedron. The important realization is that the area correction, implemented as described below, avoids this problem completely and yields extremely accurate transport properties (including the intrinsic viscosity, and for polygonal shapes), as shown for the ellipsoids in Table 1, above. Other Errors
The solution of the system of eq. (7) requires order (3N) 3 operations, and thus very long computation times could be encountered. In addition, for very large systems, the condition number for the super G matrix increases rapidly so that there may be significant round-off errors in the solution of the linear system of equations. The round-off error arises from the loss of significant figures that unavoidably occurs as the dense linear system of equations becomes large. The super matrix G is dense and nonsymmetric (even though each individual 3 ⫻ 3 block is symmetric). The round-off error, closely tied to the condition number of the matrix, effectively limits the maximum size of the matrix whose solution can be accurately computed. The presence of round-off error is detected as a nonlinearity in the graph of the transport properties vs. 1/N for N large. In single precision on a 32-bit processor, this limits the size of the matrix to about 5000 ⫻ 5000 (1700 triangles), while in double precision, one can go beyond 15,000 ⫻ 15,000 (5000 triangles). Numbers of this size are needed for the treatment of large proteins, for example ref. 39. However, single precision computations are not recommended, as discussed above. For the treatment of these large dense systems of equations, a 64-bit
*
G ⫽ r0
冤
3 ⫹ 21
冘 冘
processor is advantageous, as well as alternative techniques to perform precise computations of the interaction of elements distant from one another.40 These techniques will be implemented in an improved version of the present software. Computation of Exact Integrals of the Oseen Tensor over Triangles
In the previous section we demonstrated that it is crucial to eliminate the integration and area errors during a hydrodynamic computation to obtain high accuracy and precision. In this section we show (with details left to Appendix A), how the Oseen tensor can be integrated essentially exactly over a triangle to yield the 3 ⫻ 3 matrices of eq. (5). In so doing, we will also obtain the scaling relations for the triangle area to make the area correction. The Oseen tensor can be integrated analytically when the reference point lies within the triangle (diagonal blocks of the supermatrix), and essentially exactly when the reference point lies on a different triangle (off-diagonal elements of the supermatrix). * The results for the integrals G kj have the following form. Diagonal Blocks
The diagonal blocks are integrated exactly in the local triangle frame, and must then be rotated into the orientation they have in the laboratory frame. The local triangle frame has its z-axis along the normal to the surface, and the x-axis along one of the sides. The details are given in Appendix A, along with the definitions of the symbols shown below. The general form in the local triangle frame is,
共3dp k ⫹ 12 sc k兲
k
1 2
冘 冘 1 2
cc k
0
k
cc k
k
0
where the summations refer to the three vertex regions of the triangle. The variable r0 is the radius of the incircle of the triangle. The rotation matrix used to generate the above tensor in the laboratory frame is given explicitly in Appendix A.
3 ⫹ 21
共3dp k ⫺ 12 sc k兲
0
k
0
2 ⫹
冘 k
冥
(14)
dp k
Off-Diagonal Blocks
As detailed in Appendix A, the off-diagonal block can be represented as a nested double integral over a pair of coordinates ( x, y, varying from 0 to 1). For a triangle of area A, the general form is,
Boundary Element Method for Macromolecular Transport Properties
1199
Figure 5. Lysozyme. Space-filled model (Brookhaven file 6LYZ) and triangulation of molecular surface by MSROLL.
G ⫽ 2 A关T 0␦ ⫹ T 兴
(15)
The inner nested integral can be performed analytically with the aid of Mathematica.41 The full forms are given in Appendix A. The remaining integral in the off-diagonal case, being one dimensional, can be evaluated to any desired precision by simple Gaussian quadrature methods. To obtain integrals with 15 digits of precision, 21-point quadrature is used (the quadrature order required was determined with Mathematica, using arbitrary precision arithmetic for comparison). The integration is virtually exact in that the integration error is now much smaller than the shape, area, and discretization errors. Area Correction: Details
Previously, we pointed out that when a surface is discretized, the area decreases and the shape is deformed. The importance of the above precise relations for the integral of the Oseen tensor over a triangle is that the dependence on the area of the surface element is transparent. For the diagonal elements, the integral is effectively proportional to the square root of the area, while for the offdiagonal blocks, the integral is exactly proportional to the triangle area. It can be shown that the diagonal blocks can also be computed analytically for a spherical triangle, a computation that contains no area or shape error.42 If, on the other hand, one scales the values of the flat triangle integrals to the value of the correct area, then the scaled integral differs insignificantly from the exact integral done over a spherical triangle. This observation is the key to the area correction.
Scale the Integrals of Diagonal Blocks by f1/2 and Those of j the Off-Diagonal Blocks by fj, Where fj ⫽ Atrue /Afaceted
The algorithm thus generated essentially eliminates the integration, shape, and area errors, ill-conditioning, and allows for the very precise determination of transport properties. The Fortran program that implements the boundary element method as described above is whimsically called BEST, and it is available from the author. BEST provides a very precise full translation-rotation 6 ⫻ 6 tensor for an arbitrarily shaped, rigid molecule with stick boundary conditions. For speed of execution, the program takes advantage of LAPACK,43 which calls BLAS Level 3 machine coded routines. These are available for most processors in the market directly from the hardware vendors. The most recent release of the software also includes the computation of the intrinsic viscosity, but the implementation of that calculation will be presented elsewhere.44 Computation and Discretization of Molecular Surface Areas
How do we get the true area for a macromolecule? Connolly has developed a suite of programs to do what we need (almost!). Connolly has improved his original MS package45 by the introduction of the analytic quartic46 – 48 surface produced from the familiar process of rolling a probe sphere over the Van der Waals surface of a macromolecule. In the new MSP package,46 the surface generated is the surface that can be touched by a solvent molecule, and is called the molecular surface. The molecular surface depends on the radius of the probe sphere. If the radius of the probe sphere is zero, then the molecular surface is just the Van
1200
Aragon
•
Vol. 25, No. 9
•
Journal of Computational Chemistry
Figure 6. Linear extrapolation of the trace of the translational diffusion tensor for lysozyme vs. 1/N. N varies between 1351 and 2500 triangles. In the legend, x represents the slope of the graph, while “1” represents the intercept. The standard error is given by SE, while Tstat is the T statistic indicating the appropriateness of a linear fit.
der Waals surface produced by the intersection of the Van der Waals spheres of each atom in the molecule. Typically, the probe radius is taken to be the radius of a sphere whose area is the same as the molecular surface of a solvent molecule. For water, this value is 1.48 Å, for example. When a macromolecule can be solvated, as in the case of proteins and nucleic acids, then we must also contend with a solvation shell surrounding the molecular surface. We call that the solvent-accessible surface. The difference in the volume enclosed by the solvent-accessible surface and the molecular surface is the hydration volume. This can be computed and compared to experiment very successfully for proteins, for example. The details of that work are described elsewhere.39 Connolly’s MSP package46 and Sanner’s MSMS package49 can also produce a triangulation of the molecular surface. The raw triangulation thus produced is not suitable for boundary element calculations, however. The original intent of Connolly’s and other packages is to produce a triangulation that can be used to obtain a visually smooth representation of the surface on a computer screen. For BE calculations, one also requires that the triangles have a reasonable shape. In particular, we must avoid long slender triangles because the discretization error increases for these, and we must avoid triangles of very small area because the condition number of the supermatrix will increase. The Fortran program COALESCE has been written to address these deficiencies in either Connolly’s or Sanner’s triangulation programs, both of which produce copious amounts of small and slender triangles. This problem has been mentioned previously by Brune and Kim50 when dealing with implementations of the BE method. Allison’s implementations,31–35 on the other hand, have not attempted a
precise description of the shape of a macromolecule. Allison uses an algorithm51 which requires that all atoms or groups of atoms be represented by equal sized spheres, and thus yields only a coarse description of the molecular surface, but with very nice triangle shapes. It is also worth mentioning the SMART program of Zauhar,52 which has been constructed specifically for BE calculations and also provides very accurate surface areas. We have not used his triangulation here because of the extra complications in computing precise integrals of the Oseen tensor over his curved surface elements. It would be tempting to simplify the triangulation of the surface by limiting oneself to the Van der Waals surface instead of the solvent defined molecular surface. This approach was suggested by Zhou.53 The advantage of this method is that all surface sections are portions along a spherical surface, and this allows a comparatively simple triangulation algorithm. However, allowing solvent to enter where it physically cannot, will also yield properties such as the specific volume of a molecule that are at variance with experiment. Thus, our approach is to modify the triangulations produced by Connolly (or Sanner) to make them suitable for boundary element calculations. Our Fortran program, COALESCE, detects and processes triangles of small area, and slender acute and slender obtuse triangles. The control parameter is the fraction of the average triangle edge length that is considered to be small. This value may be used to reduce the number of triangles in a triangulation to a manageable number, at the same time that it improves the average shape of the triangles in the set. Small triangles, slender acute triangles
Boundary Element Method for Macromolecular Transport Properties
1201
Table 3. Transport Tensors for Lysozyme.
Exp. data
Calc. data
Protein
MW kDa
Hydration layer Å
D tt 10 ⫺7 cm2/s
D rr 10 7 s⫺1
D tt 10 ⫺7 cm2/s
D rr trace 107 s⫺1
Lysozyme (6LYZ)
14.3
1.1
11.2(.2)54
2.0(.1)55
11.0
2.1
(those whose two similar length edges have an acute angle between them), and slender-obtuse triangles (those whose similar length edges have an obtuse angle between them) are removed by a vertex coalescence algorithm developed by the author. This algorithm is described in Appendix B. COALESCE produces a triangulated surface of an arbitrarily shaped body with no holes (missing triangles) on the surface, taking as input, the triangulation produced by MSP or MSMS. For our application, however, Connolly’s MSP package is preferred because it computes very precise analytical surface areas required for the area correction, whereas Sanner’s MSMS does not. The triangulations with a decreasing number of triangles are easily produced from a single input data set of MSP to produce a sequence of input data for BEST and the extrapolations.
Application to Lysozyme The computations presented above for ellipsoids and benzene sufficed to illustrate the errors that needed to be addressed, and the accuracy of the computations. To demonstrate the utility of BEST and COALESCE, we apply these programs to the calculation of the transport properties of a small enzyme, lysozyme. Lysozyme is an irregularly shaped “globular” protein of molecular weight 14,300 g/mol. A triangulation of its surface is shown in Figure 5. In this case, a large number of triangles are needed to define the surface in detail, but as seen in Figure 6, the linearity of the transport coefficient graphs is excellent. The numeric precision for the determination of the intercept yields translational diffusion components to 0.02% and rotational diffusion components with numeric precision of 0.1%. In the above computations, the representation of the surface of lysozyme contains a uniform layer of hydration water 1.1 Å thick. The hydration thickness was determined by a fit to experimental data on a suite of three well-studied proteins. It is noteworthy that the high precision hydrodynamic computations will fit both the translational and rotational diffusion coefficients of the proteins within experimental error, as seen in Table 3. A thorough study of a suite of 35 different proteins, which contains a full description of this methodology, is presented elsewhere.39
Conclusions In conclusion, we have shown that a precise implementation of the boundary element method allows the computation of hydrodynamic transport tensors to a precision much higher than is acces-
sible experimentally. This tool should be able to shed light on subtle changes in transport properties that accompany conformational changes in solution due to transitions induced by temperature or salt conditions in the case of macroions, binding of substrates to enzymes, structural mutations, change of the size of the solvent, and many other effects that are buried under the typical computational errors of other widely used methods. Further extensions of the method to include the specific effects of charges and electrolytes can also be made, along the lines suggested by Allison.34 We have also presented a method, the “area correction,” which allows us to perform accurate computations even when the underlying problem is ill conditioned. This nonuniform scaling of the G matrix is receiving much greater attention in this laboratory, for therein could lie a general solution to the problem of ill conditioning in matrix inversion. In addition, the methods utilized here for the computation of exact or essentially exact integrals of the Oseen tensor over triangles can also be easily extended to other polygonal shapes and other functions. For example, if one orientationally averages the Oseen tensor, one obtains the Green function for electrostatic problems, thus the integrals for the applications of boundary element methods in those problems can also be computed essentially exactly. We have shown that an increase in the precision of the integrations pays off in greatly increased accuracy of the solutions of the dense linear set of equations. The Fortran source code, binaries, and documentation are available from the author:
[email protected].
Acknowledgments The author thanks S. Allison for helpful discussions, and T. Rosales and D. Hahn for some of the numerical computations presented in this work.
Appendix A: Integration of the Oseen Tensor over a Triangle One of the two key elements that allows the computation of the transport tensors to high numerical precision is the precise computation of the Oseen tensor over an arbitrarily oriented triangle. Let the center of this triangle be at c1 and the hydrodynamic interaction be felt at a point c2, which may be located within the triangle itself. Let the arbitrary triangle have vector vertices v1 , v2 ,
1202
Aragon
•
Vol. 25, No. 9
•
Journal of Computational Chemistry
v3. The distance that enters into the Oseen tensor [eq. (2)] is r ⫽ yi ⫺ c2, where yi is a point within the triangle. We need to use completely different strategies for the integration in the case that c2 is within the triangle itself (yielding a diagonal element of the supermatrix) and when c2 is located on a different triangle (an off-diagonal matrix element). Diagonal Elements
When the reference point c2 is within the triangle, then the distance r can be zero and the Oseen tensor will be singular. In this case, we must use polar coordinates in the local frame of the triangle to perform the integration. The natural center to take for the triangle is the incenter, and not the usual centroid. The incenter is defined as the center of the largest circle, the incircle, that will fit inside the triangle. This circle will be tangent to the triangle sides at one point each, and is shown in Figure A1. A convenient choice for the local frame of the triangle is ( xˆ ⬘, yˆ ⬘, zˆ ⬘), where the x⬘ axis is along the triangle side a, the y⬘ axis is perpendicular to it on the triangle plane, and the z⬘ axis is along the normal to the triangle. Given the vertices of the triangle, the orientation of these unit vectors in the laboratory frame is given by (vˆ 12 , nˆ ⫻ ˆ 12 , nˆ), where nˆ is the triangle normal, and ˆ 12 is the unit vector along v2 ⫺ v1, representing side a of the triangle. The triangle normal can be computed from the following cross product: nˆ ⫽ (v3 ⫺ v1 ) ⫻ (v2 ⫺ v1 )/兩(v3 ⫺ v1 ) ⫻ (v2 ⫺ v1)兩. The magnitude of the simple crossproduct gives twice the area of the triangle. The components of these vectors are the direction cosines that are needed to build the rotation matrix to compute the tensor in the laboratory frame. This rotation matrix R ⫽ [ xˆ i xˆ ⬘j ] is given by,
冋
x y z 12 12 12 共n ˆ ⫻ 兲 共n ˆ ⫻ 兲 共n ˆ ⫻ 12兲 z R⫽ 12 x 12 y nˆ x nˆ y nˆ z
册
(A1)
The tensor in the lab frame is then readily computed by: Glab ⫽ R⫺1 GlocalR. We now proceed with the integration in the local frame. The integral has contributions from two types of regions: the incircle region, and the “triangular” portions outside the incircle, the vertex regions. In polar coordinates the integrals that must be performed are the six distinct components of the symmetric local tensor (temporarily omitting the constant factor 8): G ⫽
冕冕 冋 ds r
␦ ⫹
冘冕 3
G vert ⫽
1 2
k⫽1
where
册
r r r2
ᐉk
r0
(A2)
Figure A. (1) An arbitrary triangle with an incircle of radius r 0 , an incenter c1, vertices vi . The triangle sides are a, b, c. The angle subtended by the sides a, b is C, that subtended between sides a, c is B, and that subtended by sides b, c is A. The distance from the incenter to vertex A ⫽ v3 is ᐉ a , and likewise for vertex B ⫽ v2, and C ⫽ v1. (2) A triangle with the reference point c2 outside. The vectors a and b along the sides can be used to construct a coordinate system for points inside the triangle.
We immediately note that in polar coordinates the area element ds ⫽ rdrd , thus canceling out the divergent part of the Oseen tensor. The integration over the incircle of radius r 0 is very simple and yields the tensor,
冋 册
3 G in ⫽ r0 0 0
0 3 0
0 0 2
(A3)
The vertex components can be integrated in polar coordinates as well, by collecting contributions along the arcs contained within the triangle as we vary the angle (r), whose range will vary with the distance from the incenter, r. The integration over the angle can be done immediately to yield the tensor,
冋
共1兲 共2兲 2 3共共2兲 cos2关共1兲 k ⫺ k 兲 ⫹ fk k 兴 ⫺ cos 关k 兴 共1兲 共2兲 共2兲 共1兲 2 2 3共k ⫺ k 兲 ⫺ fk dr cos 关k 兴 ⫺ cos 关k 兴 0 0
册
0 0 共1兲 2共共2兲 k ⫺ k 兲
共2兲 共1兲 共1兲 f k ⫽ 2 共sin关共2兲 k 兴cos关k 兴 ⫺ sin关k 兴cos关k 兴兲,
1
Boundary Element Method for Macromolecular Transport Properties
and (2) and (1) are the two angles that delimit the arc at a given k k value of r. The quantity ᐉ k is the distance from the incenter to the kth vertex. These quantities must be expressed in terms of the vertex coordinates. When we do so, we find that the remaining integrals can be performed exactly. Let the sides of the triangles be represented by the vectors, a ⫽ v2 ⫺ v1 ; b ⫽ v3 ⫺ v1, and c ⫽ v3 ⫺ v2. The lengths of the triangle sides a, b, and c can be computed via Pythagoras. Furthermore, let sin[C] be the sine of the angle between sides a and b, and likewise for sin[B] and sin[A]. These angles can be computed from the law of cosines. Then one can readily show that the radius of the incircle is given by,
冉
a⫹b⫺c r0 ⫽ 2
冊冉
冊
sin关C兴 1 ⫹ cos关C兴
1203
For convenience define the angles A 1 ⫽ A, A 2 ⫽ B, A 3 ⫽ C, and h k ⫽ ᐉ k /r0. The three different integrals are given by, dp k ⫽ r 0共 ⫺ A k兲共h k ⫺ 1兲 ⫺ 2共h karccos关1/hk 兴 ⫺ log关hk ⫹ sqrt关h2k ⫺ 1兴兴兲, sc 1 ⫽ ⫺
(A4)
k ⫽ 1, 2, 3
r0 cos关B ⫺ C兴兵SI1 cos关A兴 ⫹ CI1 sin关A兴其 2
sc 2 ⫽
r0 兵SI 2共1 ⫹ cos关2B兴兲 ⫹ CI2 sin关2B兴其 4
sc 3 ⫽
r0 兵SI 3共1 ⫹ cos关2C兴兲 ⫹ CI3 sin关2C兴其 4
(A8a)
(A8b)
cc 1 ⫽ ⫺r0 共sin关B ⫺ C兴cos关A兴CI1 ⫹ sin关B ⫺ C兴sin关A兴SI1 兲 and the distances of the incenter to each of the vertices are, ᐉc ⫽
cc 2 ⫽ r 0共sin关B兴2 CI2 ⫹ sin关B兴cos关B兴SI2 兲
a⫹b⫺c 2 cos关C/2兴
cc 3 ⫽ ⫺r0 共sin关C兴2 CI3 ⫹ sin关C兴cos关C兴SI3 兲
ᐉb ⫽
共a 2 ⫹ 共b ⫺ c兲 2 ⫹ 2a共c ⫺ b兲cos关C兴兲1/ 2 2 cos关C/2兴
ᐉa ⫽
共共a ⫺ c兲 2 ⫹ b 2 ⫹ 2b共c ⫺ a兲cos关C兴兲1/ 2 2 cos关C/2兴
where for k ⫽ 1, 2, 3: SI k ⫽ 2兵log关hk ⫹ sqrt关h2k ⫺ 1兴兴 ⫺ sqrt关1 ⫺ 1/h2k 兴其
sin 共1兲 ⫽ cos关␣兴, ⫽ ⫺sin关␣ ⫹ C兴,
Vertex B:
cos
共1兲
and cos 共2兲
共2兲
⫽ sin关B ⫹ ␣兴; and cos
and, 共2兲 ⫺ 共1兲 ⫽ ⫺ B ⫺ 2␣,
cos 共1兲 ⫽ sin关␣ ⫺ C兴,
and cos 共2兲 ⫽ sin关B ⫺ ␣兴,
and, 共2兲 ⫺ 共1兲 ⫽ ⫺ A ⫺ 2␣.
(A6)
With these quantities, the integrals in the tensor Gvert can now be calculated. The integrals are of three types: the angle difference, dp k , the integral over f k , sc k , and the integral over the difference of squares of cosines, cc k . The result is,
1
G vert ⫽ 2
冤
冘 k
共3dpk ⫹ 21 sck 兲
冘
cck
k
0
冘
冘
cck
0
k
共3dpk ⫺ 21 sck 兲
k
0
冘 k
0 2dpk
冥 (A7)
(A9)
It is immediately noteworthy that all components of the tensor are proportional to r 0 , or to a very good approximation, to the square root of the area of the triangle. This observation was used to generate the area correction mentioned in the text. In addition, we can note the high degree of symmetry of the formulas between vertices B and C (2 and 3) but not for vertex A, except when the integrals depend only on the difference between the two angles.
and, 共2兲 ⫺ 共1兲 ⫽ ⫺ C ⫺ 2␣,
⫽ sin关␣兴, Vertex A:
CI k ⫽ 共3 ⫺ h k ⫺ 2/h k兲
(A5)
The formulas in eq. (A4) and (A5) have the required behavior for the cases of equilateral and isosceles triangles. Define cos[ ␣ ] ⫽ r 0 /r. The angles that delimit the arcs in the vertex regions are given by, Vertex C:
(A8c)
Off Diagonal Elements
When the reference point c2 is not inside the triangle, then we can use one of the vertices as a reference point instead of the incenter (see Fig. A2). If we start at v1, then we can generate a coordinate system with coordinates ( x, y, varying between 0 and 1, and x ⫹ y ⫽ 1) that uses the vectors a ⫽ v2 ⫺ v1 and b ⫽ v3 ⫺ v1 as basis vectors. Any point inside the triangle can be found by the linear combination xa ⫹ yb, and thus, the vector to the reference point c2 from a given point in the triangle is r ⫽ v1 ⫺ c2 ⫹ xa ⫹ yb. Furthermore, the area element ds ⫽ axbdxd y ⫽ 2 Adxd y, where A is the triangle area. Thus, the tensor integral is proportional to,
冕 冕 1
G ⫽ 2 A
1⫺y
dy
0
0
dx r共 x, y兲
冋
␦ ⫹
r 共 x, y兲r 共 x, y兲 r共 x, y兲 2
册
⫽ 2 A关T 0␦ ⫹ T 兴
(A10)
It is important to note that with the given coordinate system, the integrations can be done directly in the laboratory frame and no rotation matrices are needed in this case. The integral over x can be
1204
Aragon
•
Vol. 25, No. 9
•
Journal of Computational Chemistry
performed analytically with the aid of Mathematica. Define d1 ⫽ v1 ⫺ c2 , d2 ⫽ v2 ⫺ c2, and d3 ⫽ v3 ⫺ c2, then the scalar integral is, T 0 ⫽ 兰 10 d yt 0 ( y), where,
t 0共 y兲 ⫽
冋
a 䡠 d2 ⫹ ya 䡠 c ⫹ a 冑y c ⫹ 2yc 䡠 d2 ⫹ d 1 log a a 䡠 d1 ⫹ ya 䡠 b ⫹ a 冑y2 b2 ⫹ 2yb 䡠 d1 ⫹ d12 2 2
2 2
册
The tensor integral is, T ⫽ 兰 10 d yt ( y), where, T ⫽ d 1d 1S 0 ⫹ 共d 1a ⫹ a d 1兲S 1 ⫹ 共d 1b ⫹ b d 1兲S 2 ⫹ 共a b ⫹ b a 兲S 12 ⫹ a a S 11 ⫹ b b S 22 (A12)
1
d yg 0共 y兲
S1 ⫽
0
冕
1
d yg 1共 y兲
S2 ⫽
冕
S 11 ⫽
0
冕
1
dxx 2g 2共 x兲
0
S 22 ⫽
冕
1
d yy 2g 0共 y兲
0
These integrals are performed with high precision 21-point Gaussian quadrature in BEST. The error is in the 16th decimal place! The magnitude of this error was determined by comparison to arbitrary precision computations implemented in Mathematica. Note that the process of computation of the supermatrix is an N 2 process and thus incurs an insignificant penalty compared to the N 3 process of solution of the system of equations in eq. (7). For this reason, one does not gain much in cutting corners in the integration process. The three single integrals that are required in the six integrals of eq. (A13) have the following analytical forms,
⫽
1
d yyg 0共 y兲
⫺
g 1共 y兲 ⫽
共d1 䡠 a兲2 ⫺ d12 a2 ⫺ 2yd1 䡠 ba2 ⫹ 2yd1 䡠 aa 䡠 b ⫹ y2 共a 䡠 b兲2 ⫺ y2 a2 b2 (A14)
d1 䡠 d2 ⫺ y共a 䡠 d1 ⫺ 2d1 䡠 b ⫺ a 䡠 b兲 ⫹ y2 b 䡠 c
冑d22 ⫹ 2yd2 䡠 c ⫹ y2c2
共d1 䡠 a2 ⫺ d12 a2 ⫺ 2yd1 䡠 ba2 ⫹ 2yd1 䡠 aa 䡠 b ⫹ y2 共a 䡠 b兲2 ⫺ y2 a2 b2
g 2共 x兲 ⫺
d2 䡠 a ⫹ ya 䡠 c
冑d12 ⫹ 2yd1 䡠 b ⫹ y2b2 冑d22 ⫹ 2yd2 䡠 c ⫹ y2c2
0
⫺ 冑d12 ⫹ 2yd1 䡠 b ⫹ y2 b2 ⫹
⫽
d yyg 1共 y兲
d 1 䡠 a ⫹ ya 䡠 b
0
d 1 䡠 b ⫹ xa 䡠 b
1
g 0共 y兲
and the remaining integrals are,
冕
冕
(A13)
(A11)
S0 ⫽
S 12 ⫽
d3 䡠 b ⫺ xb 䡠 c
冑d12 ⫹ 2xd1 䡠 a ⫹ x2a2 冑d32 ⫺ 2xd3 䡠 c ⫹ x2c2
共d1 䡠 b兲2 ⫺ b2 d12 ⫺ 2xd1 䡠 ab2 ⫹ 2xd1 䡠 ba 䡠 b ⫹ x2 共a 䡠 b兲2 ⫺ x2 a2 b2 (A16)
Note that these integrals are given in terms of vector inner products, and can be easily computed from the coordinates of the triangle vertices.
Appendix B: The Algorithms in COALESCE The basic input into COALESCE is a set of vertices and other information provided by a triangulation program such as MSP or MSMS. COALESCE first computes the lengths of all the triangle sides in the set, and finds the average length. The program then searches for all triangles whose sides are smaller than the length of interest. This length is given as a fraction of the average triangle side upon input. Thus, f ⫽ 0.1 implies that all triangles whose sides are less than 10% of the average are declared small and will be coalesced out of existence. When COALESCE finds such a pair of triangles (a side is always shared with another), it proceeds to
(A15)
superimpose the vertices across from the small side, coalescing these vertices, and the small triangles are eliminated from the set. This is shown schematically in Figure B1 for the case of a slender-acute triangle. Vertices a and b of the small triangles are superimposed. This implies that all the triangles with vertex at a will have their new vertex assigned as b. The vertex with the least number of triangles sharing it is the one that is moved. Vertex b is shared with four other triangles (other than the small ones), while vertex a is shared with only three triangles. This algorithm will remove small triangles regardless of their shape, and will remove triangles that are slender-acute when only one side is small. Slender-obtuse triangles are handled differently. These triangles do not have any sides that are small. These are reshaped by an algorithm that leaves the number of triangles in the set the same, as shown in Figure B2. COALESCE also checks for and eliminates other rare instances such as dangling triangles (those who share only one side with other members of the set), and triangles of zero area (degenerate triangles). Because the process of eliminating slender-acute triangles can produce other slender-acute triangles, the processing of a triangulation is done in two passes. The program produces two output files, a vertex file to serve as input to BEST, and a histo-
Boundary Element Method for Macromolecular Transport Properties
Figure B. (1) Slender-acute triangle coalescence. (2) A slender-obtuse triangle at the bottom is eliminated by partitioning the large triangle above it from vertex a to vertex b.
gram file with information about the processing and the triangle size distribution. The algorithm in COALESCE scales as N 2 .
References 1. Giacovazzo, C., Ed. Fundamentals of Crystallography; Oxford University Press: New York, 1992. 2. Richards, E. G. An Introduction to Physical Properties of Large Molecules in Solution; Cambridge University Press: New York, 1980. 3. Berne, B.; Pecora, R. Dynamic Light Scattering: With Applications to Chemistry, Biology and Physics; Wiley–Interscience: New York, 1976. 4. Eden, D.; Elias, J. G. In Measurement of Suspended Particles by Quasi-Elastic Light Scattering; Dahneke, B., Ed.; Wiley–Interscience: New York, 1983. 5. Stryer, L. Science 1968, 162, 526. 6. Swaminathan, R.; Hoang, C. P.; Verman, A. S. Biophys J 1997, 72, 1900. 7. Ryba, N. J. P.; Marsh, D. Biochemistry 1992, 31, 7511. 8. Sanders, J. K. M.; Hunter, B. K. Modern NMR Spectroscopy; Oxford University Press: New York, 1987. 9. Harding, S. E.; Rowe, A. J.; Shaw, W. V. Biochem Soc Trans 1987, 15, 513. 10. Ediger, M. D. Annu Rev Phys Chem 1991, 42, 225. 11. Aragon, S. R.; Pecora, R. Macromolecules 1985, 18, 1868. 12. Happel, J.; Brenner, H. Low Reynold’s Number Hydrodynamics; Prentice Hall: Englewood Cliffs, NJ, 1965.
1205
13. Kim, S.; Karilla, S. J. Microhydrodynamics; Butterworth-Heinemann: New York, 1991. 14. Bloomfield, V. A.; Dalton, W. O.; van Holde, K. E. Biopolymers 1967, 5, 135; Ibid. Biopolymers 1967, 5, 149. 15. Garcia de la Torre, J.; Bloomfield, V. A. Q Rev Biophys 1981, 14, 81. 16. Garcia de la Torre, J.; Bloomfield, V. A. Biopolymers 1977, 16, 1779. 17. Teller, D. C.; Swanson, E.; de Haen, C. Methods Enzymol 1979, 61, 103. 18. Oseen, C. W. Hydrodynamik; Academiches Verlag: Leipzig, 1927. 19. Yamakawa, H. J Chem Phys 1970, 53, 436. 20. Rotne, J.; Prager, S. J Chem Phys 1969, 50, 4831. 21. Wilson, R. W.; Bloomfield, V. A. Biopolymers 1979, 18, 1205. 22. Garcia Bernal, J. M.; Garcia de la Torre, J. Biopolymers 1981, 20, 129. 23. Antosiewicz, J.; Porschke, D. J Phys Chem 1989, 93, 5301. 24. Zipper, P. Analytical Ultracentrifuge Users Meeting, University of Nottingham, 1998. 25. Pastor, R. W.; Karplus, M. J Phys Chem 1988, 92, 2336. 26. Venable, R. M.; Pastor, R. W. Biopolymers 1988, 27, 1001. 27. Garcia de la Torre, J.; Huertas, M. L.; Carrasco, B. Biophys J 2000, 78, 719. 28. Youngren, G. K.; Acrivos, A. J Fluid Mech 1975, 69, 377. 29. Hu, C. M.; Zwanzig, R. J Chem Phys 1974, 60, 4354. 30. Wegener, W. A. Biopolymers 1986, 25, 627. 31. Allison, S. A.; Nambi, P. Macromolecules 1992, 25, 3971. 32. Allison, S. A.; Nambi, P. Macromolecules 1994, 27, 1413. 33. Allison, S. A.; Tran, V. T. Biophys J 1995, 68, 2261. 34. Allison, S. A.; Mazur, S. Biopolymers 1998, 46, 359. 35. Allison, S. A. Macromolecules 1999, 32, 5304. 36. Ladyzhenskaya, O. A. The Mathematical Theory of Viscous Incompressible Flow; Gordon & Breach: New York, 1963. 37. Kirkwood, J. G.; Riseman, J. J Chem Phys 1948, 16, 565. 38. Brenner, H. J Colloid Interface Sci 1967, 23, 407. 39. Aragon, S.; Rosales, T.; Hahn, D. K. Biophys J, in preparation. 40. Zauhar, R. J.; Varnek, A. J Comp Chem 1996, 17, 864. 41. Mathematica, a computer algebra package by Wolfram Research, Champaign, IL. 42. Aragon, S., unpublished work. 43. Anderson, E.; Bai, Z.; Bischof, C.; Blackford, S.; Demmel, J.; Dongarra, J.; Du Croz, J.; Greenbaum, A.; Hammarling, S.; McKenney, A.; Sorensen, D. LAPACK User’s Guide; SIAM: Philadelphia, 1999, 3rd ed. 44. Hahn, D. K.; Aragon, S. R., in preparation. 45. Connolly, M. L. QCPE Bull 1981, 1, 75. 46. Connolly, M. L. J Mol Graphics 1993, 11, 139. 47. Connolly, M. L. J Appl Crystallogr 1983, 16, 548. 48. Connolly, M. L. Science 1983, 221, 709. 49. Sanner, M. F.; Spehner, J. C.; Olson, A. J. Biopolymers 1996, 38, 305. 50. Brune, D.; Kim, S. Proc Natl Acad Sci USA 1993, 90, 3835. 51. Pakdel, P.; Kim, S. J Rehol 1991, 35, 797. 52. Zauhar, R. J. J Comput Aided Mol Des 1995, 9, 149. 53. Zhou, H. Biophys J 1993, 65, 955. 54. Sophianopoulos, A. J.; Rhosdes, C. K.; Holcomb, D. N.; van Holde, K. E. J Biol Chem 1962, 237, 1107.
