INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS APPLICATIONS CONCLUDING REMARKS COMPUTER PROGRAMS REFERENCES

Many properties of dilute solution of macromolecules, or dilute suspensions of colloidal particles, are related to the Brownian motion of rigid particles constituting the solute. In a dilute system (as it will be assumed throughout this paper), the diffusivity of the particles is solely determined by their size and shape, along with the temperature and the solvent viscosity. The Brownian trajectories followed by the rigid particles have a random nature, but the statistics of the properties associated to those trajectories is determined by a 6 × 6 diffusion tensor, D, which, as described in the next section, contains terms related to the translational and rotational diffusivities (Brenner, 1967; Garcia Bernal and García de la Torre, 1980; Harvey and García de la Torre, 1980). Fundamental theories describe how the translational displacement (Einstein, 1906) and the rotational motions (Favro, 1960) of the Brownian particle can be predicted in terms of D. In contrast, D can be calculated from the shape of the rigid particle. Years ago, such calculation was only possible for simple, symmetric shapes such as revolution ellipsoids (Perrin, 1934, 1936; Koenig, 1975) or cylinders (Tirado and García de la Torre, 1979, 1980). Nowadays, methodologies that enable the accurate calculation of D of arbitrarily shaped particles have emerged. A possible approach is the boundary element method (Youngreen and Acrivos, 1975a,b; Allison and Nambi, 1992). We also have the technique of bead modeling (Kirkwood, 1954; Bloomfield et al., 1967; García de la Torre and Bloomfield, 1981) (for a recent review, see Carrasco and Garcia de la Torre, 1999), which is implemented in useful, publicly available computer programs (Carrasco and Garcia de la Torre, 1999; Garcia de la Torre et al., 1994, 2000a,b). Thus, from all the existing knowledge, the Brownian diffusivity of the isolated particle, which does not experience interactions with other particles or external fields, can be adequately described and predicted.

The problem is more complex in some situations of great relevance, in which the particles in the dilute solution interact with some external agent. A number of experiments could be mentioned, including polar or polarizable macromolecules in external fields (Vesterberg, 1993; Chae and Lenhoff, 1995; Hagerman, 1996; Allison, 2001), encounter-complex formation (Northrup et al., 1984; Gabdoulline and Wade, 1997), or particles diffusing in the presence of impenetrable (and perhaps electrified) obstacles or macroparticles (Huertas et al., 1996, Fernandes et al., 2001). In these situations, the Brownian trajectories and the observable properties depend, in an intricate manner, not only on the single-particle D tensor, but also on the specific features of the environment. Even for the simplest cases, theoretical results will not be available, and the proper approach has to be a computer simulation.

In this work, we describe a Brownian dynamics simulation algorithm, implemented in the computer program BROWNRIG, for rigid particles of arbitrary shape that interact with an external agent. In the design of the computer program, a separate module, specific of each system, evaluates the forces and torques experienced by the particle under the action of an external agent. Other modules, common to all cases, use those forces and torques along with the diffusion tensors to evaluate the Brownian displacements. Thus, the procedure can be adapted, with a minimum effort, to cases with arbitrarily complex environments.

There are some antecedents of our work. Allison (1991) described a method with body-fixed and lab-fixed frames, but his work was limited to free diffusion and mainly intended for polarized light scattering. Some authors (see for instance, Wade, 1996) have devised simulations in which the diffusional anisotropy of the particles is ignored; this is tantamount to preaverage orientationally the diffusion tensor, which is thus like that of an equivalent sphere. This approach may be particularly inadequate for nonglobular, aspherical particles. In our procedure, the particle is regarded as an arbitrarily shaped body, with fully anisotropic translational and rotational diffusivities, and translation-rotation coupling.

In the second part of this paper, we present several applications of this methodology, including free diffusion, transport in an electric field, and motion in the vicinity of barriers. These examples serve not only to test the performance of BROWNRIG, but they also illustrate the potential applications of the procedure and describe how to analyze the simulated trajectories to calculate the observable properties.

	THEORY AND METHODS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS APPLICATIONS CONCLUDING REMARKS COMPUTER PROGRAMS REFERENCES

Generalities

The evolution of the macromolecular particle during its Brownian trajectory is characterized by some vectorial or matricial quantities that change with time. For a vectorial quantity v, the time dependency will be indicated as v(t). The Cartesian components of v(t) a given time t will, in turn, depend on the system of reference axes. The ultimate frame or reference will be a system of laboratory-fixed axes (LFAs), but we also need a system of particle-fixed axes (PFAs), which, obviously, will change as the particle moves. When necessary, the axis to which the components of v(t) is referred will be denoted by a superscript. So, v^(l)(t) indicates components of v(t) referred to the LFAs, whereas v^(t₀)(t) indicates the components of v at time t, referred to PFAs that the particle had at some other time t₀. The same notations will be used for matricial quantities. It is obvious that, if v is an intrinsic property of the particle (e.g., the diffusion tensor, the dipole moment, the coordinates of some of its points, etc.), then its components at any time referred to the PFAs at that time will be always the same: v^(t)(t) = v^(t₀)(t₀) = v⁽⁰⁾(0) or, simply, = v.

It is convenient to place the particle in the initial instant (t = 0) of the simulation with its center of diffusion, r, at the origin of the LFAs so that r(0) = 0. We recall that the center of diffusion is the proper point to which translational and rotational quantities must be referred (García Bernal and García de la Torre, 1980; Harvey and García de la Torre, 1980; García de la Torre and Bloomfield, 1981). Its position is one of the results of the usual bead-model calculations (García de la Torre et al., 1994). Also, at this initial instant, the PFAs coincide with the LFAs; in other words, the zero-time particle axes coincide with the lab axes, and therefore subscripts (l) and (0) mean the same.

The transformation of the coordinates of a vectorial quantity or the components of a matricial quantity in a system of axes (a) into those for another system of axes (b) is made through a transformation matrix M^(ab),

(1)

where (a) and (b) can be any of the PFAs or the LFAs. These matrices are orthogonal, and therefore M^(ab) = M^(ba)T. The superscript T indicates matrix transposition or column vector. In the formalism of our Brownian dynamics algorithm, we have to handle matrices that transform from the PFAs at some time t to the LFAs, M^(tl).

The output of the simulation program, i.e., the trajectory, will consist of successive registers, each containing the coordinates of the center of diffusion at time t referred to the LFAs, r^(l)(t), and the matrix M^(tl), which transforms from the PFAs at time t to the LFAs. If v is some specific particle-fixed vector, its coordinates in the LFAs at time t will be

(2)

Eq. 2 is written taking into account that, at any time, v^(t)(t)  v. If v is one of the unitary vectors of the PFAs, e_k, then its coordinates at time t in the LFAs system, e(t), are the elements of the kth column of M^(tl).

Brownian dynamics

For the purpose of numerical simulation, the Brownian trajectory is discretized as a series of time steps of identical duration, t. At each step of its trajectory, the particle experiences a Brownian translational-rotational displacement that can be described by a generalized vector of dimension 6:

(3)

where x(t₀), k = 1, 2, 3 (x, y, z) are translational displacements starting at time t₀, referred to the particle coordinate frame at t₀. _k, k = 1, 2, 3 are rotations around the t₀ particle's frame axes, e_k(t₀).

The displacement has two components:

(4)

The component

(5)

is the Brownian displacement. It obeys the generalized Einstein's law,

(6)

where D_lm are the components of the 6 × 6 generalized diffusion tensor, D, with 3 × 3 boxes D_tt, D_tr = D and D_rr (Brenner, 1967; Harvey and García de la Torre, 1980). Here, subscripts t and r mean translation and rotation, respectively:

(7)

If a force F or a torque T are acting on the particle, they make a contribution to its displacement,

(8)

The mobility tensor that determines this displacement is (k_BT)¹D, where k_B is the Boltzmann constant and T the absolute temperature. If we define a generalized force, , with

(9)

where we specify the initial instant at which they apply, and refer them to the PFAs of the initial instant, then

(10)

The generalized 6 × 6 diffusion tensor D will be the object of a separate, previous calculation. As commented upon in the introduction, for a rigid particle of arbitrary shape, this quantity can be obtained from rigid-body hydrodynamic calculations, based either on the boundary element methods (Youngreen and Acrivos, 1975a,b; Allison, 1999) developed by other authors, or on bead- or shell-model procedures developed in our group and implemented in public-domain programs like HYDRO (García de la Torre et al., 1994) and its extensions (García de la Torre et al., 2000a,b; 2001; Garcia de la Torre, 2001a; Garcia de la Torre and Carrasco, 2002).

The computational simulation will also require the evaluation, at the beginning of every time step, of the forces, F^(t₀)(t₀), and torques, T^(t₀)(t₀), acting on the particle. The corresponding calculation will be specific of the field or external agent that the particle is experiencing. A frequent case will be that of polar or polarizable particles in an electric field. The joint description of the electrostatics and hydrodynamics of a particle bearing charges, surrounded by its ionic atmosphere, may be remarkably complex. At any rate, our BROWNRIG program is devised modularly, so that any procedure for F^(t₀)(t₀) and T^(t₀)(t₀) can be implemented as a separate subprogram, and linked with the rest of the software.

As an oversimplified example of electrostatic interaction, valid in some instances and useful as an illustration of the procedures, is the case in which the particle has n_q charges q_s, s = 1, ... , n_q, rigidly attached at some points within the particle. The particle is in a position-dependent electric field E(r), whose value at the point occupied by charge s at time t₀ is E(t₀) = E^(l)(r(t₀)), with

(11)

where r is the position vector of the center of mass, and r_sc is the vector going from the center of mass to the point where charge s is located. Then, the force acting on charge s is F(t₀) = q_sE(t₀), and the total force, referred to the PFAs at t₀ is

(12)

The torques at each point charge are referred to the center of diffusion, so that

(13)

Algorithm and trajectory simulation

A previous phase in the calculation is the evaluation, using HYDRO or related software (García de la Torre, 1994; 2000a), of the diffusion tensor with reference to the PFAs, D (= D⁽⁰⁾(0) = D^(t₀)(t₀) for any t₀). The HYDRO calculation also provides the location of the diffusion center, which will be the origin of the PFAs at any time, and the point to which origin-dependent quantities are referred. Initially, at t = 0, we obviously have

(14)

and

(15)

where I is the 3 × 3 identity matrix.

At a given time, t₀, from which a Brownian simulation step is taken, the position vector of the center of the mass and the matrix of transformation to the LFAs are r^(l)(t₀) and M^(t₀l), respectively. The displacements q^B (Eq. 5) and q^d (Eq. 10) are added (Eq. 4) to obtain the generalized displacement vector, q, or more properly denoted q^(t₀)(t₀). Its first three components give the translational displacement of the center of mass, r^(t₀)(t₀) = (x(t₀), x(t₀), x(t₀)). This displacement is added to the previous position of the center of mass to obtain the new one; for this purpose, r^(t₀)(t₀) must be previously transformed to the LFAs, and thus we have

(16)

The rotational step consists of three successive rotations of the particle around the axes 1, 2, and 3 (fixed as they stand at time t₀). The orientation of any particle-fixed vector v changes from v(t₀) to v(t) after the step. The new coordinates in the t₀ axes are given by

(17)

where R₁, R₂, R₃ are the matrixes that transform coordinates upon axial rotation, for instance,

(18)

where c₁ = cos ₁ and s₁ = sin ₁, and they have the property R = R, k = 1, 2, 3. From the definition of the M^(tl) matrixes, we can write that, v^(l)(t) = M^(tl) · v^(t)(t), and also v^(l)(t) = M^(t₀l) · v^(t₀)(t). Combining these expressions with Eq. 17, we obtain

(19)

Eqs. 16 and 19 are the basis of the Brownian dynamics algorithm, because they give the coordinates of the center of mass, and the transformation matrix into the LFAs, of the particle-fixed axes at time t = t₀ + t from their values at time t₀, with t being the time step.

The three components and the three angles ₁, ₂, ₃ of r^(t₀)(t₀) are random numbers with zero mean and a covariance matrix equal to 2Dt, according to Eq. 6. They can be generated as described, for instance, by Allison and McCammon (1984).

	APPLICATIONS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS APPLICATIONS CONCLUDING REMARKS COMPUTER PROGRAMS REFERENCES

Free diffusion

In experimental situations, the translational diffusivity will be measured by the translational diffusion coefficient, D_t. According to Einstein (1906), this quantity determines the time-dependence of the mean-square displacement of the center of mass of the particle as

(20)

where r_c(t₀) is the position vector of the center of mass, with respect to a laboratory-fixed system, at time t₀, and r_c(t₀ + t) is the position when some time, t, has elapsed. The average r_c(t₀) · r_c(t₀ + t)_t0 is made over all the possible choices of the initial instant, t₀, along the Brownian trajectory. The theory of translational diffusivity states that the translational diffusion coefficient is given by D_t = Tr(D_tt)/3, where Tr(D_tt) is the sum of the diagonal components of tensor D_tt (Eq. 7).

Rotational diffusion is experimentally observed from the time (or frequency) dependence of some dynamic electro-optical or spectroscopic properties, measured with techniques such as fluorescence anisotropy (Belford et al., 1972), electric birefringence and dichroism (Wegener et al., 1979) and NMR relaxation (Woessner, 1962) (see García de la Torre et al. (1997, 1999, 2000b), connecting these properties with rigid-body hydrodynamics). These techniques monitor the Brownian reorientation of certain vectors within the particle; more specifically, the observed properties are related to the time dependence of the function, briefly denoted as P₂(t), defined by

(21)

Actually, P₂(t) is the average (over the choice of initial state) of the second Legendre polynomial of the angle subtended by two successive orientations of some particle-fixed unitary vector, u(t₀) and u(t₀ + t), separated by time t. Eq. 21 is valid for Brownian motion in general circumstances, even in the presence of fields or restricting effects. If we consider the case of free diffusion, then further results are available that allow a full description of the rotational dynamics. For free diffusion, P₂(t) is a multiexponential decay function having up to five components,

(22)

where the five relaxation times, _k, are essential, dynamic properties of the particles, solely determined by the components of the rotational tensor, D_rr (Eq. 7). If D₁, D₂, and D₃ are the eigenvalues of D_rr, the relaxation times are (in ascending order) ₁ = (6D_r 2)¹, ₂ = (3(D_r + D₁))¹, ₃ = (3(D_r + D₂))¹, ₄ = (3(D_r + D₃))¹, ₅ = (6D_r + 2)¹, where D_r = (D₁ + D₂ + D₃)/3 and  = (D + D + D  D₁D₂  D₂D₃  D₁D₃)^1/2. The five amplitudes, a_k, depend on the components of tensor D_rr and on the coordinates of vector u in the system of eigenvectors of D_rr; the expressions are simple but lengthy and can be found elsewhere (see, for instance, García de la Torre et al., 1999). In summary, for free diffusion, the rotational dynamics can be entirely predicted using Eq. 16 from the D_rr tensor, which is in turn calculated by rigid-body hydrodynamics.

To test the performance of the present simulation methodology, the trajectories generated by BROWNRIG can be analyzed and compared with the results predicted by the rigid-body hydrodynamic theory. The molecule that we analyzed was the human antibody IgG3. This molecule, whose bead-model, depicted in Fig. 1, was proposed by Gregory et al. (1987), has three regions: two Fab and one Fc joined together by a long hinge. We have chosen this molecule because, beside its intrinsic interest, it is the example model that comes along with the HYDRO package (http://leonardo.fcu.um.es/macromol) and can be implemented quite straightforwardly. As can be seen in Fig. 1, the IgG3 model consists of 15 beads, and this is the model we have used to calculate the hydrodynamic properties of IgG3. The output of this calculation, performed with HYDRO, namely the coordinates of the center of diffusion and the generalized diffusion tensor, is used as the input of BROWNRIG. Doing so, we simulated a 1 × 10⁹-step trajectory of IgG3, with a time step of 1 × 10¹³ s, at 298 K in a solvent with viscosity of 0.01 Poise.

View larger version (125K): [in this window] [in a new window]   FIGURE 1   Bead model of the human antibody IgG3, showing the two Fab, the Fc, and the hinge regions.

The translational diffusivity of IgG3 can be characterized by the mean square displacement of the coordinates of its center of mass, and is described by Eq. 20. Analyzing this property, we found that the (isotropic) translational diffusion coefficient of IgG3 recovered from our simulation is 2.293 × 10⁶ cm²s¹, which is in excellent agreement with the value of 2.288 × 10⁶ cm²s¹ for the same coefficient computed by HYDRO.

We also studied the rotational dynamics of IgG3 computing P₂(t) from the simulated trajectory according to Eq. 21. The P₂(t) function is related to the decay of fluorescence anisotropy, r(t), and also with the nuclear magnetic resonance relaxation. The unitary vectors used (u(t₀) and u(t₀ + t) in Eq. 21) were the ones that coincide with the direction of the segments that join the beads modeling the Fab and Fc regions. The results are presented in Fig. 2, where P₂(t) functions resulting from simulation are depicted and compared to P₂(t) calculated theoretically using Eq. 22 (the value of D_rr was calculated with HYDRO). It can be seen that the agreement is excellent at least until the P₂(t) decays to 0.1, and, after that, the results are accompanied with noise but still agree within simulation error.

View larger version (17K): [in this window] [in a new window]   FIGURE 2   P2 function decay of the unitary vectors that superimpose the segments that join the beads modeling the Fab and Fc regions of IgG3.

Limited diffusion

An example of limited diffusion concerns the situation of a particle embedded in a lipidic bilayer. The model system we have used to simulate this situation consists of an ellipsoidal particle, namely a prolate ellipsoid with long and short semi-axes b = 7.5 Å and a = 5 Å, respectively, with axial ratio b/a = 1.5. For a prolate ellipsoid, the 6 × 6 generalized diffusion tensor, D, is diagonal with components D = 3.97 × 10⁶ cm²s¹, D = 3.67 × 10⁶ cm²s¹ (double) and D = 1.00 × 10⁹ s¹, D = 7.36 × 10⁸ s¹ (double) (Perrin, 1936; Koening, 1975). The movements of the ellipsoid are restricted along the z axis by two infinite planes (here viewed as the membrane interfaces) separated by distance h as depicted in Fig. 3. This way, and depending on the position of the center of mass of the ellipsoid, some orientations are allowed and others are not, because the ellipsoid cannot trespass the membrane interfaces (if, in a step, the ellipsoid trespasses the walls, then the ellipsoid is placed in its previous position and the step is repeated with a new random displacement). An equilibrium property that can be defined to characterize the orientational restriction imposed by the two planes is the order parameter, S, defined as
where  is the angle formed by the long molecular axis and the laboratory-fixed z axis, and P₂(cos )  (3 cos²  1)/2 is the second Legendre polynomial of cos . The order parameter S is obtained by calculating P₂(cos ) at each step and averaging over the entire trajectory. S can assume values between 0.5 (for a molecular axis orientation perpendicular with respect to reference axis z) and 1 (for a molecular axis orientation parallel with respect to reference axis z). A particle with no orientational restrictions will have a value for S of 0. P₂(t), as defined in Eq. 21, is related with S if we notice that P₂(0) = 1 and, when t  , P₂(t)  P₂². Because S = P₂, we have that P₂(t) decays from 1 to S². We can see in Fig. 4 how the P₂(t) of the unitary vector along the long axis of the ellipsoid decays from 1 until it reaches a baseline of 0.03. These results were obtained for a separation of planes, h, of 15 Å. The value of S for the same ellipsoid and in the same conditions, calculated as an average over the entire trajectory, is 0.181 (S² = 0.03 and is in agreement with the value of the baseline of P₂(t). This numerical result for S is in good agreement with the result (S = 0.201) given by an analytical expression derived for the orientation of ellipsoids induced by planar obstacles (Fernandes et al., 2001). The decay profile of P₂(t) and the value of S will depend on the geometrical parameters of the simulated system. As h becomes bigger compared to the dimensions of the particle (which means less orientational restrictions), S will tend to 0 and will approximate the value of particle experiencing free diffusion (results not shown). As discussed in our previous paper (Huertas et al., 1996), the short time rotational dynamics is characterized by an initial relaxation time given by 1/_ini = [d lnP₂(t)/dt]_t=0, which must coincide with that of a freely diffusing particle, _free. From the numerical simulation results, as can be seen in Fig. 4, we obtain the initial slope _ini = 4.52 × 10⁹ s¹, in excellent agreement with the value for free diffusion of the ellipsoid _free = 6D = 4.41 × 10⁹ s¹ given for the dynamics of ellipsoids (Perrin, 1936; Koening, 1975).

View larger version (13K): [in this window] [in a new window]   FIGURE 3   Representation of the system used to simulate limited diffusion: z is the distance from the center of the ellipsoid to the plane, a and b are the values for the short and long semi-axes, respectively, and h is the separation between the planes that limit the diffusion.

View larger version (9K): [in this window] [in a new window]   FIGURE 4   P2(t) decay of the unitary vector along the long axis of the ellipsoid (b = 7.5 Å; a = 5 Å). Results were obtained for a separation of planes, h, of 15 Å. The straight line describes the short time rotational dynamics of the ellipsoid.

The translational properties of the ellipsoid are also characterized by the mean square displacement of the coordinates of the center of mass. Regarding Eq. 20, that gives the general law of the translational diffusion, we will have to distinguish two situations. The diffusion along directions x or y, which are not bounded, obeys the law presented in Eq. 20 (the constant 6D_t is substituted by 2D_t because we are considering diffusion along a single direction). However the diffusion along the z direction, which is bounded, does not follow Eq. 20. The diffusivity of the ellipsoid along the z direction is presented in Fig. 5. It can be seen that, for short times, the ellipsoid will not have, on average, its diffusion limited by the planes and presents a linear behavior according to Eq. 20 (Huertas et al., 1996). For longer times, the ellipsoid will have its diffusion, along the z direction, limited by the planes, and the mean square displacement will reach a plateau (Huertas et al., 1996) as can be seen in Fig. 5. At long times, z(0) and z(t) are not correlated and [z(t)  z(0)]² = 2z²  2z². The asymptotic value in the correlation is 3.10 Å², whereas, for the z² and z² values obtained as averages over the entire trajectory, we reached 2.98 Å², which is in good agreement.

View larger version (11K): [in this window] [in a new window]   FIGURE 5   Translational dynamics of the ellipsoid (b = 7.5 Å; a = 5 Å; planes separated by h = 15 Å) given by the square displacement of its mass center.

Application of an external electric field

Let us now consider the dynamics of a charged, rigid macromolecule in the presence of an external electric field. To simulate this situation we first considered a system consisting of a rigid molecule of ten identical spherical subunits, having a diameter of 20 Å, forming two arms connected by a joint, with the arms making an angle of 30 degrees (see Fig. 6 A). The model molecule is submitted to a constant uniform electric field along the z direction. This electric field has a constant positive value everywhere in the space. The internal charges are placed in the extremity subunits of each arm and at the tangency point of the two central units.

View larger version (10K): [in this window] [in a new window]   FIGURE 6   Structure of the model molecule, a rigid bent rod, used to simulate the dynamics with an applied external electric field. The black dots represent the electric charges. (A) The angle between the arms is 30°. (B) The angle between the arms is 120°.

We simulated the electrophoretic migration of charged particles in an applied external electric field using the above-described molecule with a central charge of +3e and two extremity charges of 1e. Thus, the molecule has a net positive charge +1e, and a dipole moment directed along the bisector of the interarm angle (axis y). In Fig. 7 A, we present the results of the electrophoretic migration, given by z coordinates of the diffusion center, along the direction of the applied electric field. We first consider an applied field of 1.00 × 10⁶ N · C¹, which is sufficiently weak so that the Brownian tumbling of the particle overcomes the orientational effect arising from the interaction of the electric field with the dipole moment, and the particle migrates with no preferential orientation, a situation similar to what happens in experimental systems. The progression of the z coordinate with time can be compared with the expected theoretical displacement, represented as a straight line in Fig. 7 A. In general, the particle moves with a velocity given by v = M · E, where E is the applied electric field and M is the electrophoretic mobility tensor, which depends on both hydrodynamic and electrostatic properties of the particle and its environment. For our simple representation of electrophoretic migration, the velocity is given by v = Eq_net/f_eff, E is the applied electric field, q_net is the net charge of the molecule and f_eff is the effective friction coefficient, which is related with the components of the friction tensor  and depends on the degree of orientation of the particle. Because, in this situation, the particle does not migrate with any particular orientation, f_eff is simply kT/D_t and can be obtained from the output of HYDRO. Substituting the values of electric field, net charge, and effective friction coefficient, we reach a value of 2.715 × 10³ m · s¹ for the migration velocity. Figure 7 A shows that the agreement between the simulated trajectory and that expectation is rather good; the fluctuations around the straight line are simply a manifestation of the Brownian drift that is superimposed to deterministic migration.

View larger version (16K): [in this window] [in a new window]   FIGURE 7   Electrophoretic migration of a charged particle, given by the coordinates of its mass center, submitted to an applied external electric field of (A) 1 × 106 N · C1 and (B) 5.5 × 107 N · C1.

In Fig. 7 B, we present the results of the electrophoretic migration of the same particle, now for a much stronger applied field of 5.5 × 10⁷ N · m¹. It can be seen that the coordinates show an almost linear progression versus time, and the simulation results fit very well the theoretical result given by the straight line. We used the same expression to calculate the theoretical migration velocity as in the previous situation. The difference is that, in this case, the applied electric field is very intense and the particle is strongly oriented when migrating, with the particle-fixed y axis nearly aligned with t, the lab-fixed z axis. In this situation, f_eff is the _yy of the translational friction tensor, which, from the output of HYDRO, is found to be f_eff = 5.36 10¹¹ Kg · s¹. Substituting the values in the expression, we reach a value of 0.165 m · s¹ for the migration velocity. Two points are worth noticing: 1) the applied electric fields are intense, particularly in the second case. We did so, knowing that it had no correspondence with a real situation, to obtain a clear electrophoretic migration without consuming a lot of computation time; and 2) the representation of the electric properties of the molecule is very simple, and we neglected the simulation of effects such as the counter-ion atmosphere (Allison, 1992; Allison, 2001). We did so because we just wanted a straightforward quantitative illustration of the general performance of our algorithm.

To simulate the orientational dynamics of a dipole in an applied external electric field, manifested in techniques such as transient electric birefringence or dichroism, we used the model molecule depicted in Fig. 6 B now with 11 subunits and the two arms making an angle of 120°. The molecule has a permanent dipolar moment: the central unit has a charge of +2e and the two extremity units have a charge of 1e each, which means a zero net charge. The applied electric field has a magnitude of 5.5 × 10⁷ N · C¹. In Fig. 8, we present the results for the time-dependent normalized electric birefringence, defined as
where n_sat would be the birefringence reached in an infinitely strong field, P₂ is the second Legendre polynomial, ₁ and ₂ are the angles that the first and second arms of the molecule make with the direction of the applied electric field, respectively. Using the Brownian dynamics algorithm, the individual trajectories of a thousand molecules were simulated. At any given time, P₂(cos ₁) and P₂(cos ₂) are calculated for each molecule, and, from the whole sample average, we extract the value of the birefringence. The time course of the birefringence is illustrated in Fig. 8. When the electric field is applied, the birefringence goes down until it reaches a plateau value, n' = 0.10, which agrees well with the prediction from theoretical expressions for the steady-state electric birefringence of bent-rods (Mellado and Garcia de la Torre, 1982). When the field is switched off, the birefringence returns to zero. The jumps from n = 0 to the steady value n' = 0.10 and back to n = 0 are visualized simply in Fig. 8. A detailed analysis of the kinetics of the decay and rise processes (although not included within this work) is obviously feasible. Therefore, our methodology seems to be a promising tool for the simulation of transient electric birefringence or dichroism of rigid particles with arbitrary shapes in electric fields of arbitrary intensity.

View larger version (16K): [in this window] [in a new window]   FIGURE 8   Profile of the average electric birefringence of an ensemble of 1000 molecules with a time varying-applied external field.

	CONCLUDING REMARKS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS APPLICATIONS CONCLUDING REMARKS COMPUTER PROGRAMS REFERENCES

We have presented a Brownian dynamics algorithm that allows the calculation of the trajectory of the anisotropic diffusion of an arbitrary shaped rigid particle that interacts with other particles or external fields. This algorithm uses the general diffusion tensor obtained from HYDRO (Garcia de la Torre et al., 1994) to evaluate the displacements experienced by the diffusing particle. This way, we can obtain a more realistic dynamics and, consequently, more accurate results can be withdrawn from simulations. This algorithm can be used in a variety of situations (like the example applications we have presented and many others) where the arbitrary shape of a particle, interacting with external agents, is determinant for its dynamics and properties.

	COMPUTER PROGRAMS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS APPLICATIONS CONCLUDING REMARKS COMPUTER PROGRAMS REFERENCES

We will make the source code of the routine that implements BROWNRIG algorithm publicly available from our web site (http://leonardo.fcu.um.es/macromol/). The eventual user has to design a personal module to evaluate the forces and torques experienced by the particle under the action of an external agent. With this modular procedure, the user can describe the interactions of the particle with the desired degree of sophistication needed for the user's particular system and situation.