Stochastic Dynamics of Bionanosystems: Multiscale Analysis and Specialized Ensembles
Stephen Pankavich, Yinglong Miao, Jamil Ortoleva, Zeina Shreif, Peter Ortoleva
SStochastic Dynamics of Bionanosystems: Multiscale Analysis and Specialized Ensembles
S. Pankavich a , Y. Miao b , J. Ortoleva, Z. Shreif b , and P. Ortoleva b a Department of Mathematics b Department of Chemistry Center for Cell and Virus Theory Indiana University Bloomington, IN 47405
Abstract
An approach for simulating bionanosystems such as viruses and ribosomes is presented. This calibration-free approach is based on an all-atom description for bionanosystems, a universal interatomic force field, and a multiscale perspective. The supramillion-atom nature of these bionanosystems prohibits the use of a direct molecular dynamics approach for phenomena like viral structural transitions or self-assembly that develop over milliseconds or longer. A key element of these multiscale systems is the cross-talk between, and consequent strong coupling of, processes over many scales in space and time. Thus, overall nanoscale features of these systems control the relative probability of atomistic fluctuations, while the latter mediate the average forces and diffusion coefficients that induce the dynamics of these nanoscale features. This feedback loop is overlooked in typical coarse-grained methods. We elucidate the role of interscale cross-talk and overcome bionanosystem simulation difficulties with (1) automated construction of order parameters (OPs) describing supra-nanometer scale structural features; (2) construction of OP dependent ensembles describing the statistical properties of atomistic variables that ultimately contribute to the entropies driving the dynamics of the OPs; and (3) the derivation of a rigorous equation for the stochastic dynamics of the OPs. As the OPs capture hydrodynamic modes in the host medium, “long-time tails” in the correlation functions yielding the generalized diffusion coefficients do not emerge. Since the atomic scale features of the system are treated statistically, several ensembles are constructed that reflect various experimental conditions. Attention is paid to the proper use of the Gibbs-hypothesized equivalence of long-time and ensemble averages to accommodate the varying experimental conditions. The theory provides a basis for a practical, quantitative bionanosystem modeling approach that preserves the cross-talk between the atomic and nanoscale features. A method for integrating information from nanotechnical experimental data in the derivation of equations of stochastic OP dynamics is also introduced.
Keywords : virus, ribosome, bionanosystems, time-of-flight experiments, chemical tagging, nanopore devices, multiscale analysis, coarse-graining, Gibbs hypothesis, Smoluchowski equation, long-time tails 1
I. Introduction
Viruses, ribosomes, and other bionanosystems (BNS) are supramillion atom structures that must be understood in terms of the interplay of atomistic, highly fluctuating behaviors and more coherent dynamics involving the collective motion of many atoms simultaneously. The challenge we address is to quantify this picture and use it to understand BNS behaviors such as structural transitions, self-assembly, and interaction with other features (e.g., membranes) in a cell’s interior or other biological milieu. Being supramillion atom in size, a BNS presents a grand challenge for predictive modeling. For example, the highly optimized molecular dynamics code NAMD, run on a 1024 processor supercomputer, can simulate approximately 1.1 nanoseconds of evolution for a complete satellite tobacco mosaic virus (STMV) in one day. However, viral structural transitions take a millisecond or longer. Thus, using this approach would take about 2500 years. While some other methods taking advantage of the symmetry of the BNS (e.g., icosahedral viruses) can proceed much faster, they cannot address the highly nonlinear, local and dissipative nature of the BNS phenomena of interest. Phenomenological models (e.g., wherein peptides and nucleotides are represented as beads, or other major subunits of the BNS are lumped into simplified computational elements) must be recalibrated with each new application, severely limiting their predictive power. Lumped models can also lead to difficulties in that the definition of the subunits may not be invariant over the time course of the phenomena of interest. Thus, new concepts and computational algorithms are needed for understanding and designing bionanosystems. Advances in multiscale theory imply an algorithm that enables all-atom dynamical computer simulations of a BNS over biologically relevant time periods. Here we present a rigorous and generalized reformulation of the all-atom multiscale approach. Multiscale analysis is a way to study systems that simultaneously involve processes on widely separated time and length scales. It has been of interest at least since the work on Brownian motion by Einstein . In these studies, FP (Fokker-Plank) and Smoluchowski equations are derived either from the Liouville equation or via phenomenological arguments for nanoparticles without internal atomic-scale structure. We extended this work by (1) accounting for atomic-scale internal structure of the BNS; (2) introducing general sets of structural order parameters (OPs) characterizing nanoscale features of the system; (3) inventing a way to introduce OPs into the analysis without the need for tedious bookkeeping to ensure that the number of degrees of freedom is unchanged; and (4) introducing ensembles constrained to fixed values of the OPs to construct average forces and diffusion coefficients in the coarse-grained equations of stochastic dynamics . These elements underlying our BNS modeling approach are strongly interrelated.
Order Parameters
Major features of a BNS change slowly in time relative to the 10 -14 second characteristic timescale of fast atomic collisions and vibrations. Variables describing these nanoscale features include the CM (center-of-mass), orientation, and nanometer-scale BNS substructural units (e.g., capsomers for a virus). OPs change slowly for several reasons: (1) they have large inertia as they describe the coherent motion of many atoms simultaneously; (2) fluctuating atomic forces tend to cancel as they are averaged over the BNS surface; (3) frictional forces are large at the nanoscale and thereby repress rapid motions; and (4) energy barriers may be large relative to thermal energies. As a result, although OPs have some macroscopic character, they represent features which are small enough that many associated phenomena must be understood in terms of a stochastic description. 2
Ensembles and Interscale Cross-Talk
While the OPs change slowly, the atomistic variables explore a broad range of configurations. As one is usually interested in the dynamics of the OPs, and has little concern for the detailed atomic configurations at any instant, it is relevant and practical to treat these atomistic variables using a probabilistic framework. To do so, we construct a probability density that accounts for the instantaneous values of the OPs at each point along their time course. Like total system energy for the canonical ensemble used in Gibbsian statistical mechanics of macroscopic systems, values of the OPs must be used to construct ensembles characterizing the likelihood that the system resides in each configuration of the atomistic variables. As for the thermodynamics of macroscopic systems, the construction of the ensembles developed in Section II reflect the available information on system preparation and the exchange of energy and mass during the process of interest. Unlike macroscopic systems, however, one may not always choose the ensemble that is most convenient for carrying out the calculations, that is, the finite size of the system requires greater attention when constructing the ensemble and choosing among several possible ensembles, as noted in Sect. II.
Perturbation Theory and Coarse-Grained Stochastic Order Parameter Dynamics
One may identify small parameters (notably mass, time, length ratios, and force strengths) that can be used to develop an orderly perturbation solution of the Liouville equation for the N -atom probability density. The lowest order solution is constructed via an OP-constrained entropy maximization principle and a closely related restatement of the Gibbs hypothesized equivalence of ensemble and long-time averages. Expansion methods based on the smallness of these ratios imply equations for the reduced probability density describing the stochastic evolution of the OPs. Depending on the character of the system and the physical/biological regime of interest, the stochastic equation derived is of the Smoluchowski or Fokker-Planck type. The frictional and average force factors appearing in these equations depend sensitively on the values of the OPs and on the conditions (e.g., iso-energetic or isothermal) to which the system is subjected. Our multiscale perspective provides insight into the cross-talk among the atomistic and nanoscale variables. As suggested in Fig. 1 , the OPs set the overall (nanoscale) context that determines the relative probability of atomistic configurations, which in turn, mediate the average forces and frictional effects that drive the dynamics of the OPs. Such a relationship is the basis of the feedback loop expressed in
Fig. 1 . This is in sharp contrast to typical coarse-grained models, which aim to simplify the system by reducing the number of interacting elements. This is usually done via “decoupled” coarse-graining wherein atoms are lumped together, and an effective force on the lumped element is set forth . In particular, a systematic method to derive the parameters needed for the effective force field from trajectories and force data collected from a single MD simulation was developed . This approach was applied to various systems such as liquids , monosaccharides , and peptides . The results were in good agreement with atomistic simulations while computational cost was significantly reduced. Note however, that MD simulations can only be carried over short time and length scales. Thus, one wonders if the derived parameters would still be applicable over longer scales. Furthermore, the derived parameters are used as a basis of the coarse-grained simulation without regard for the feedback loop of Fig. 1 (i.e., the role of the overall structure on the computation of the coarse-grained forces is overlooked). For example, the ensemble of fluctuations of atoms deep within a tightly wound globular protein is likely to differ significantly from that when the protein is uncoiled. A similar comment is appropriate for the transition state between two distinct conformations separated by a high energy barrier. A further difficulty is that this type of coarse- 3 graining provides no self-consistency criterion regarding the number of atoms to be used in a lumped element. The ensemble of atomic-scale fluctuations that should be used in a proper computation of the coarse-grained forces depends on temperature and other conditions, thus the parameters in the coarse-graining should change with these conditions. Trajectories produced in this manner are only reliable over short times wherein the overall nanostructure doesn’t change enough to significantly alter the ensemble of atomistic fluctuations. The feedback loop of
Fig. 1 suggests the need for the simulation algorithm in
Fig. 2 . The coarse-grained equations for stochastic OP dynamics are solved numerically. As the average forces and diffusion coefficients in these equations change with the OPs, they must be co-evolved with the OPs. These factors are expressed as statistical averages weighted with OP-constrained probability densities.
Fig. 2
In the fully coupled multiscale simulator, the stochastic equations of the Langevin type for the BNS structural OP dynamics are simulated. Diffusion coefficients and average forces are computed on-the-fly as they change with the OPs via the feedback loop of Fig. 1.
Input BNS structure BNS structural time-course Construct diffusion coefficients Stochastic simulator Set biophysical conditions and run parameters Construct average forces Output data analysis results
Average forces and diffusion coefficients
Order parameters
Ensemble of atomistic configurations
Fig. 1
OPs characterizing nanoscale features affect the relative probability of the atomistic configurations which, in turn mediates the forces driving OP dynamics. This feedback loop is central to a complete multiscale understanding of nanosystems and the true nature of their dynamics.
II. Ensembles for Atomistic Fluctuations of a BNS
A BNS is envisioned to have dual macroscopic/microscopic characteristics - nanoscale components, each consisting of many atoms, and atomic-scale, highly fluctuating features. To characterize the coherent nature of the nanoscale features, OPs are introduced and expressions relating them to the all-atom configuration are developed (see Section III and appendix). These OPs transcend particular atomic associations; for example, the CM of a nanoparticle is a collective property of thousands of atoms, not of a single atom. Involving so many atomic notions in a coherent way, OPs in a BNS typically change slowly in time relative to the timescale of individual atomic vibrations and collisions (e.g., 10 -14 to 10 -12 seconds). Thus, associated with the coherent nanometer-scale, stochastic atomic-scale dichotomy, there is a corresponding timescale separation. The objective of the formulation developed here is to simultaneously account for these diverse scales and the cross-talk among them. A common theme of many statistical mechanical theories of nanosystems is that the timescale separation suggests that the atomistic variables explore a representative set of configurations over a period of time on which the OPs are relatively constant. The relative residence time in each member of this ensemble of configurations is characterized by a probability density . By definition, quantifies the likelihood that the short-scale behaviors visit each configuration in this ensemble. Given our lack of knowledge of most of the atomic-scale detail, and that most of it is not directly relevant to understand BNS behavior, stability, and function, an entropy-maximization approach is used to construct . We develop this approach and integrate it with multiscale analysis in this and subsequent sections. In the statistical mechanics of macroscopic systems, this enables the study of various experimental conditions (e.g., isothermal versus iso-energetic). Here, we show that this can be used to introduce nanotechnical measurement information into the analysis. Given knowledge of some nanoscale features of the system derived from nanotechnical measurements, our objective is to construct the probability distribution for the rapidly fluctuating degrees of freedom. The statistical state of a BNS for most biologically relevant conditions evolves slowly, i.e., with a characteristic time much greater than that (10 -12 to 10 -14 seconds) of atomic vibrations and collisions. In Section IV, we show that this implies that the probability density characterizing this statistical state must depend on the (the set of 6 N atomic positions and momenta) only via the set of M structural OPs plus the energy . In contrast, if depends on in other ways, then evolves on the atomic time scale. Thus, while the lowest order analysis of the N -atom system implies that has the restricted form , it does not provide further information on the specific functional form of . Since no information about the detailed atomic-scale state of the system is usually known, we turn to an information theoretic method to construct . 5 According to Jaynes , the entropy S for the classical system as formulated here for the nanosystem is given in terms of an integral of . For the restricted state counting as in Appendix A, this becomes ( II.1 ) where is a product of
M+1
Dirac delta functions, one for each of the M structural OPs and one for the energy. The construction of follows from a consideration of the particular experiment under investigation. Several distinct cases are analyzed below.
Iso-energetic Closed Systems
Consider an isolated system in which the energy is held constant at a value E . To construct , we maximize S subject to normalization . ( II.2 ) Since energy is automatically fixed to via the -function mediated state counting, no further constraints on S maximization are required here. Using the Lagrange multiplier approach, one obtains the probability density denoted in the present case by in the form ( II.3 ) . (
II.4 ) Using (
II.1 ) and (
II.2 ) one finds that . Since gradients of Z with respect to the OPs will be shown in Section IV to derive coherent OP dynamics, one concludes that in this iso-OP ensemble, evolution is driven by entropy differences. This is an expression of the Second Law for an isolated nanosystem. We consider to be a conditioned probability density. If is the probability that the system is in a state where the M +1 OPs are in a small element about the particular value (i.e., ), then the corresponding probability density is . Here and in the following, we drop the for simplicity as all the equations for (as in Section IV) are linear so that such factors are found to cancel in the equations. With the above, the overall probability for this case is given by ( II.5 ) a general form that, with modification of the factors, holds for all cases considered in this study. A superficial examination of (
II.4 ) suggests that Z is not bounded. However, integration over the 3 N atomic positions is bounded by the volume of the vessel containing the system. Furthermore, the delta function on energy restricts the momentum integrations because the total energy is a constant, the potential energy is bounded from below (i.e., the potential energy can never approach as the core atom-atom potential becomes as two atoms overlap), and the energy is a monotonically increasing function of all the atomic momenta. Isothermal Closed Systems
Consider an ensemble of systems in contact with a thermal bath so that the energy is only known by its average , , (
II.6 ) 6 where is the set of M structural OPs, is a particular set of their values, and (see Appendix A). With this average energy constraint and that of normalization, entropy maximization implies ( II.7 ) . (
II.8 ) The isothermal partition function Q is related to the free energy via . This result was obtained in an earlier study in a more circuitous manner. Mixed Energy Systems
In two stage experimental protocol, the system is initially equilibrated with a thermal bath, and allowed to evolve in an energy-conserving vessel. For example, a nanoparticle first traverses a thermalizing gas and is then redirected to a vaccine chamber for iso-energetic evolution. The statistical analysis for such an experiment proceeds in two related stages. First, one solves the problem iso-energetically arriving at for particular energy value E . Taking the energy reference state such that the statistical average of any quantity in this mixed ensemble takes the form ( II.9 ) i.e., the problem is solved for each particular energy E and the result is then averaged via the canonical weight and . Nanotechnical Experiments
The ensembles constructed via the entropy maximization procedure as above are designed to experimentally integrate derived information with the theory. For example, the classic canonical ensemble imposes the temperature to construct the ensemble of atomistic configurations while the microcanonical ensemble imposes the energy to arrive at the ensemble. In a similar manner, we treat experimental approaches wherein the apparatus is designed to constrain the system so that the uncertainty in the spectrum of atomistic fluctuations is decreased. As the allowed spectrum of atomistic configurations and fluctuations determines the average forces and diffusion coefficients determining the stochastic OP dynamics, this experimental information can greatly improve our capability to predict BNS dynamics.
1. Fluorescence Tagged Ensembles
Fluorescence or mass labeling can be used to identify parts of a BNS which reside near its surface during the process of interest . The label is presented via the microenvironment and only penetrates within a short distance below the BNS surface. The notion of a molecular surface has previously been quantitatively defined . This definition provides an algorithm for constructing a surface denoted here such that if the atomic configuration is , then points for which lie within the structure and those for which are in the microenvironment. In our tagged ensemble approach, we first construct the average depth below the surface where the atoms in the labeled parts of the BNS reside. For arbitrary , will be on the order of magnitude of the size of the BNS. By experimental design, it is known that is small for configurations during the process of interest. Thus, to eliminate irrelevant 7 configurations, we limit the count of states by adding a factor to all integrals over , where is a depth to which the labeling compound penetrates below the BNS surface. The tagging of near-surface components of a BNS yields information about the system; thus, we expect to have a related decrease in entropy S , . ( II.10 ) The additional information has the effect of decreasing the entropy and thereby altering the driving force for OP dynamics. By construction, S is an increasing function of . Thus, if the tagging is only very close to the surface and these tagged BNS subunits are always near the surface, then any configuration with these subunits far below the BNS surface would have a low probability and thereby be avoided by the stochastic OP dynamics we derive (see Section IV). Though the BNS remains with the tagged subunits near the surface in one experiment, they may not remain so in another. Thus, a simulation that showed a divergence of the tagged entropy descending to large negative values because of energetic driving forces would indicate a structural inversion (i.e., inside-out) transition.
2. The TOF Ensemble
Time-of-flight experiments measure the cross-section for collision between the background gas molecules and the BNS. Efficient algorithms exist for computing , the cross-section for the specific microconfiguration of the BNS . Thus, if is the observed value, then to the counting factor we add , which is unity when is within of and zero otherwise. For example, when the OPs cause D to approach zero, the TOF entropy approaches . As the average force on the OPs is the gradient of the TOF entropy with respect to the OPs, this implies that the dynamics of the TOF ensemble member systems be driven to evolve away from such states. As shown in Section IV, the Langevin equations we develop will guarantee that if the OPs are initially consistent with , they will remain so for the simulated time-course.
3. The Nanopore Ensemble
Confined BNS evolution in a nanopore is being studied for several reasons : (1) the BNS is fixed in space so that it can be subjected to electrical, dielectric, or other fields and mechanical disturbances, and (2) the size of the pore can be chosen to only allow the entry of nanoparticles in a given size range. Furthermore, once the BNS is totally fitted into a nanopore, zero electric forces can be induced in the pore material and forces thereby applied to induce structural transitions. The presence of the nanopore is already accounted for in our formulation as it is reflected in the limits on the positions of the atoms – i.e., they remain in the system which is within the nanopore. A goal of multiscale analysis is to construct the coarse-grained probability density . For an experiment represented by one obtains ( II.11 ) This states that W is the reduced probability density for the N -atom distribution as constructed here. However, as the multiscale analysis develops, we find that the form does not capture corrections due to the slow evolution of the OPs described in the Liouville equation. 8 In our earlier treatment , we also constructed , but in two steps. We present that approach in Appendix B to contrast that presentation to the present one. To proceed with the multiscale analysis, we develop OPs suited for a spectrum of BNS (Section III) and then use them with a multiscale analysis of the Liouville equation to determine the time dependence of the statistical state of these systems (Section IV). The properties of the ensembles constructed affect the character of the stochastic OP dynamics. One such property is the average of the momentum of a single atom ( i here). The probability density in all the above cases depends on via a term in the kinetic energy. This term is symmetric (even) in while itself is asymmetric (odd). Thus, involves the integral over a symmetric range in each component of (i.e., for the x -component of ) with for the isothermal case and for an isolated system at energy E . This implies . Thus, the average of any linear combination of the is zero, a conclusion shown in Section IV to affect the form of the stochastic equation for OP dynamics. III. OPs for Bionanosystems
A key element of multiscale analysis of a BNS is the identification of OPs that describe its nanoscale features. The multiscale analysis requires that these parameters must have several characteristics. They should be flexible, meaning, applicable to a system of multiple nanoscale components such as viruses, functionalized nanoparticles, or cell membranes. They must capture the hydrodynamic modes of the host fluid to avoid “long-time tail” difficulties in the construction of diffusion coefficients that arise in the theory of Brownian motion. The set of OPs must be “complete”, i.e., they do not couple to other slowly evolving variables omitted from the set. For practical reasons, OPs should also be implementable in an automated computational procedure. In contrast, one might divide the system into “structural units” (e.g., protein complexes or viral capsomers). Unfortunately, such structural units are often artificial, i.e., are more geometric than a natural consequence of the interatomic forces and the level of thermal agitation. Finally, a central characteristic of an OP is that it evolves slowly relative to the timescale of atomic vibrations and collision. In a sense these parameters define the multiscale character of the N -atom system. In Section IV we show that failure to base a multiscale analysis on such OPs is at the heart of the weakness of other coarse-grained approaches. The first studies to observe such behaviors were conducted by both Rahman and by Alder and Wainwright using molecular dynamics. Rahman determined the velocity autocorrelation function in liquid argon, while Alder and Wainwright showed that the velocity autocorrelation function of a Brownian particle in a hard-disc and hard-sphere fluid decays as , where is the dimensionality. Long time tails were explained by a hydrodynamic model wherein the movement of the Brownian particle creates a vortex in the host fluid which in turn affects the behavior of the Brownian particle. This constitutes an additional delayed effective interaction between the particle and its host medium. Hence, hydrodynamic modes in the host medium create memory effects and long-time tails in the correlation function which can cause anomalies in the diffusion coefficients. Since the work of Alder and Wainright, there have been several attempts to derive the asymptotic time behavior of the velocity autocorrelation function . This was also shown 9 experimentally using light-scattering , diffusive wave spectroscopy , and neutron scattering . For historical record, however, it is important to mention that the first hydrodynamic theory of translational Brownian motion was provided by Vladimirsky and Terletzky in 1945. These authors derived an equation that is equivalent to that used later for the mean square displacement of the Brownian particle . To start our multiscale analysis and to include OPs that capture hydrodynamics and other slow modes, we set forth a methodology based on a deformation of space that our recent studies suggest is ideally suited even for complex BNS phenomena . A central property of these or other OPs is that they evolve slowly. Slow OP dynamics emerge in several ways including: • by inertia associated with the coherent dynamics of many atoms evolving simultaneously; • in the case of migration over long distances; • via stochastic forces that tend to cancel; • in species population levels as in chemical kinetics or self-assembly which involve many units, only a few of which change in an interval of time relative to the collision time. These factors are realized for a number of variables: • structural parameters characterizing objects such as viral capsids or ribosomes that stay intact during a transition; • orientational angles for a nanostructure; • density-like variables characterizing the hydrodynamic modes or composition; • scaled positional variables describing the motion of disconnected molecules across a nanostructure such as an oil droplet; • curvilinear/twist parameters describing the conformation of RNA, DNA and proteins . While the OPs and factors that make them evolve slowly require approaches differing in technical details, the methods presented here can be extended to them all. To illustrate our methodology, we focus on intact nanostructures. OPs for this case are introduced by embedding the system in a volume . Basis functions for a triplet of labeling indices are introduced which are orthonormalized. If computations are carried out using periodic boundary conditions to simulate a large system (e.g., to minimize boundary effects and to handle Coulomb forces), periodic basis functions can be used. Here, the nanosystem deforms in 3-D space. Points in this space are considered to be a displacement of original points . As presented in our earlier study , a set of vector OPs are constructed using orthonormal polynomials in atomic coordinates via the procedure of Appendix C: ( III.1) where is the mass of atom ( ) and the integral orthogonality of the basis functions implies that the matrix is nearly diagonal. With these OPs, a multiscale Molecular Dynamics/Order Parameter extrapolation (MD/OPX) approach has been developed to simulate large BNS . In the approach, a short MD run estimates the rate of change of the OPs, which is then used to extrapolate the system over times that are much longer than the 10 -14 second timescale of fast atomic vibrations and collisions. It has been shown that the OPs satisfy all criteria for a multiscale computational approach and the simulator has provided one to two orders of magnitude computational speed-up over direct MD. 10 IV. Multiscale Theory of BNS Kinetics
An equation of stochastic OP dynamics is now obtained via a multiscale analysis for a classical N -atom system that preserves the feedback between the atomistic and nanoscale variables of Fig. 1 . The rapidly fluctuating degrees of freedom are hypothesized to explore a representative sample of configurations during a fraction of one characteristic time for the dynamics of the OPs. Using the Gibbs hypothesis, an average over such a time interval is equated to an average over the representative sample of configurations of the rapidly fluctuating variables. In Section II such an ensemble is constructed to be consistent with given values of the slowly changing OPs. Thus, such ensembles are a key element of our multiscale BNS theory. In this section we use these concepts and the Liouville equation to derive equations for the slow stochastic dynamics of a BNS. The analysis starts by writing the Liouville equation for the N -atom probability density , i.e., for Liouville operator . Many authors (see [3, 4] for reviews) have recast this equation in the form for small parameter (e.g., a mass ratio as in Section III). This recast equation is solved perturbatively via a Taylor expansion in . Recently , we have argued that the above multiscale form of the Liouville equation arises by making the ansatz that depends on the N -atom state both directly, and indirectly via a set of OPs, . In adopting this perspective, is not a set of additional independent dynamical variables, rather its appearance in is a place-holder for a special dependence of on that underlie the slow temporal dynamics of . This perspective avoids the need for tedious bookkeeping wherein a number of atomistic variables are removed to preserve the total number ( N ) of degrees of freedom. The dual dependence of on can be ensured if is sufficiently small. In turn, this is assured if is slowly varying in time relative to the timescale (10 -14 to 10 -12 seconds) of atomistic fluctuations. Because many BNS processes depend sensitively on their internal atomistic structure, and earlier studies on nanoparticles via the Liouville equation ignored this internal structure, we develop an approach that is symmetric with respect to all atoms in the system (e.g., those in a virus, cell membrane, and aqueous medium). In the following, we derive an equation for the OP probability distribution in two cases: (1) an isolated system of given energy and (2) an isothermal system. The mixed and other cases of Section II are addressed briefly. A. The Isolated System at Energy E
The coarse-grained probability density for a system of energy E with structural OPs is defined via . (IV.1 ) The operators and noted above arise out of the ansatz on the dual dependence of on (i.e., ) and the chain rule. Thus, involves partial derivatives with respect to at constant (when operating on in the multiscale form ), and conversely for . Again, this does not imply that the are dynamical independent variables. Rather, this is a reflection of the dual dependence of However, as the system is isolated, 11 has no derivative with respect to . This occurs as the walls of the vessel are perfectly reflecting, that is, conserving of kinetic energy and redirecting velocity. As shown in the analysis of the OPs (Section III), the parameter (<<1) naturally emerges from Netwon’s equations, i.e., when computing for momentum-like variable . Here we use it to organize the construction of . In Section III, was found to be the ratio of small-to-large characteristic masses (e.g., the mass of a typical atom to that of the major nanoscale components of a BNS). For simplicity of notation, we label the OPs such that are the structural OPs of Section III and indicates the energy H . As with other multiscale approaches reviewed in Section I, it is postulated that depends on the sequence of times where . The times for are introduced to account for the slower behaviors in , while accounts for processes on the fast timescale (i.e., changes by one unit when 10 -14 seconds elapse). With the above framework, the Liouville operators take the form ( IV.2 ) (
IV.3 ) where is the set of all OPs except energy H . There is no term in since the conjugate momentum for H is zero due to the iso-energetic nature of the walls continuing the system. Note, and are the mass, momentum, position and net force for atom i . Since neither involves a derivative with respect to H , only depends on energy parenthetically. If E is the fixed value of H over the evolution in the iso-energetic ensemble of interest here, then the probability density is denoted . The strategy we now pursue is to construct via an asymptotic expansion in and use it to construct an equation for the evolution of . The Gibbs hypothesized equivalence of long-time and ensemble averages plays a key role in our multiscale analysis of the Liouville equation. The lowest order Liouville operator has a corresponding propagator that evolves dynamical variables (notably any function of ) in time, yet leaves unchanged. Thus, the Gibbs hypothesis takes the form . ( IV.4 ) Here indicates an average of any variable over all configurations selected by the factor and weighted by the conditional probability for , i.e., (
IV.5 ) This provides key elements needed to carry out a multiscale analysis for iso-energetic systems. We seek a perturbative solution of the Liouville equation in the form with quasi-equilibrium character to lowest order, i.e., is independent of and arises form entropy 12 maximization. To lowest order in we find . (IV.6)
Recalling results from Section II for the iso-energetic case and that any function of is in the nullspace of , one obtains (
IV.7 ) for and factor to be determined. To the multiscale Liouville equation implies . (
IV.8 ) Taking to be at , this equation has the solution . (
IV.9 ) Since the system is bounded in space by the walls of the vessel, and as the potential energy approaches when any two atoms overlap, for any function of and fluctuates but remains finite for all . With this, there is no term to balance the divergent contribution to and hence must vanish. A general equation for the coarse-grained probability density for any of the ensembles of Section II can be obtained by using the Liouville equation. By definition, is related to the N -atom probability density via ( IV.1 ). The Liouville equation implies (IV.10) for N -atom density Properties of the delta function ∆ , and integration by parts imply . (IV.11) Developing in a series in , one can construct the ( n +1)-th correction to the rate of change of from the n -th order correction to . For the OPs of Section III and the triple index labeling used there, one obtains . (IV.12) As noted in Section II, the ensemble average of any of the individual atomic momenta is zero, therefore so is the iso-energetic ensemble average of . The factor in , the lowest order solution for , is seen to be the lowest order contribution to in . With the above, we see that the contribution to is zero, and hence must be . Using the expression for and the definition of , it is found that the only contribution to is from . Notice that , the initial data for , is arbitrary, 13 i.e., depends on the experiment of interest. Thus, for some initial data, it is seen that the contribution to can have short timescale dependence (e.g., due to a shock wave). However, for most BNS phenomena we do not expect such phenomena to be part of the experimental design. This implies that for the BNS phenomena of interest, is in the null space of , i.e., depends on only through . The result is that in order for (IV.11) to be closed in to , must be a function of only. For the special case , up to , and this is the special case studied henceforth, i.e., we assume the system is initially in the quasi-equilibrium state. Collecting the above results for the special case with iso-energetic conditions implies a Smoluchowski equation for : (
IV.13 ) (
IV.14 ) (
IV.15 ) As and (see Section II), it is seen that for the closed iso-energetic system. Hence, OP dynamics is entropy-driven. This is a natural result for an isolated system. The above analysis can be extended to higher order in . If all initial data appearing in higher-order corrections to (i.e., the ) are zero, closure of the equation for is ensured. This does not create difficulties in resolving secular behavior in the higher order analysis as discussed in detail elsewhere . Corresponding to the Smoluchowski equation derived above there is an equivalent set of friction-dominant (i.e., non-inertial) Langevin equations for the stochastic dynamics of . These equations provide a convenient way to simulate a BNS via an algorithm as in
Fig. 2 . As and depend on they must be computed at each timestep. Thus, the feedback in
Fig. 1 is accounted for.
B. Isothermal Systems
The development for a system at temperature closely follows that of Section IV.A except that the lowest order solution is given by . (
IV.16 ) It must be recognized that this is a heuristic treatment. For example, a system in contact with a thermal bath experiences continuous exchange of energy through the bounding walls so that in a 14 more rigorous treatment the Liouville equation must be solved with stochastic boundary conditions. Again, the objective is to derive an equation for the probability density, which for this case is given by (
IV.17 ) for isothermal N -atom probability density . One proceeds in a similar manner as above but with subtle differences. Using the expression for , we find for free energy see Section II.B. As before, we find a Smoluchowski equation for up to : , ( IV.18 ) , (
IV.19 ) . (
IV.20 ) This result is similar to that for the iso-energetic case except that the diffusion coefficients and average forces are modified to reflect the ensemble of atomistic fluctuations for the versus the ensemble. Unlike for the thermodynamic limit , these factors could be quite different for the two ensembles.
C. Initial Isothermal, Iso-energetic Evolution Cases
Consider the two stage experiment wherein the system is first equilibrated with a thermal bath and subsequently allowed to evolve iso-energetically. Since the Liouville equation is linear, a superposition of solutions also satisfies the equation. Thus, the solution corresponding to a linear combination of states with different initial probability densities, each of which is thereafter iso-energetic, can be summed with a weighting by the probability that the system had a given energy initially. As noted in Section II, the probability density for the energy of systems equilibrated with a bath at temperature is (for the energy conversion ). Therefore (
IV.21 ) where subscript mix indicates the two stage experiment. This result shows that are related via Laplace transformation. However, the average forces and diffusion coefficients for the iso-energetic case cannot simply be averaged via the weight and used in a Smoluchowski equation for . Rather, these factors must be used in the Smoluchowski equation for and the solution transformed to obtain .
D. Nanotechnical Experiments
Ensembles were considered in Subsection II.D for various nanoscale experiments. They were (1) chemical labeling to distinguish subunits of a nanostructure on its outer surface; (2) 15 confinement of a BNS in a nanopore; and (3) selecting nanoparticles of a given cross-section area throughout the process of interest. Each case has a distinct partition function and thereby average force and diffusion coefficients that mediate the evolution of the OPs. The Langevin equations for the three cases can be used to study structural transitions of a BNS under the given conditions. Ensembles of the various types could reveal particular aspects of the system. For example, a labeling experiment could reveal results on swelling transitions in a viral capsid since the penetration depth of the labeling reflects the diffusibility of the label molecule in the swollen versus normal state. Inclusion of step function factors in the state counting of Section II allows for a situation wherein there is a range of values for which the partition function is zero. This implies that as this region of –space is approached, will diverge, driving the system away from this region. This is how the step function condition in the expression for the partition function guides the Langevin evolution away from configurations prohibited by the experimental design (e.g., a virus cannot swell beyond the volume of a nanopore in which it is trapped.
V. Conclusions
A methodology to obtain Smoluchowski equations for the stochastic dynamics of OPs characterizing key nanoscale features of a BNS is presented. Unlike decoupled coarse-graining methods, the key feedback of
Fig. 1 accounting for the modification of the average forces on the instantaneous values of OPs is accounted for in our fully coupled multiscaling strategy. Our development begins with the automated construction of OPs that are demonstrated via Newton’s equations to evolve on timescales long relative to that of atomistic fluctuations. Probability densities are developed to characterize the statistics of rapidly fluctuating atomistic degrees of freedom. The densities explore a representative set of configurations during a fraction of the characteristic time on which the OPs evolve. The atomistic degrees of freedom and the OPs are captured by the dependence of the aforementioned probability density on the OPs and, conversely, the contribution of the rapidly fluctuating to the entropy and the determination of OP dynamics by averaged energy. In this sense, our methodology can be termed fully coupled multiscaling. This is in contrast to decoupled coarse-graining where the effect of the evolving OPs on the computation of the forces on lumped elements is ignored. In our approach, the Liouville equation is solved perturbatively upon identification of a smallness parameter , that is a ratio of characteristic masses, lengths, times, or energies. The result is a Smoluchowski equation for the stochastic dynamics of the OPs. The average force and diffusion coefficients that appear in this equation depend on the OPs. Contrastingly, in a decoupled method, coarse-grained forces are computed without regard to the OP-dependent averaging that is always changing as, for example, a viral structural transition unfolds. A reconsideration of ensembles of atomistic fluctuations (Section II) allows one to tailor the Gibbs hypothesis to a variety of special cases that introduce modern nanotechnical experimental approaches. These include fluorescent or mass labeling, TOF selected nanosystems, and nanopore confined experiments. Our use of the Gibbs hypothesis, rather than integration over atomistic degrees of freedom as in other multiscale approaches, enables full theory/experiment integration. As the systems of interest are nanometer in scale, the usual equivalence of results of different ensembles (e.g., the microcanonical and canonical) may not hold; rather, these 16 differences may reflect distinct experimental protocols, and alter the resulting average forces and diffusion coefficients that appear in the Smoluchowski equation. Our multiscale analysis suggests that such a Smoluchowski equation can only be derived when the perturbation parameter is sufficiently small. Thus, coarse-graining is not meaningful unless the lumped variables for which forces are computed evolve significantly slower than the atomistic ones. Thus, coarse-grained potentials for amino acids or nucleotides may not provide a viable starting point as they are still rapidly fluctuating variables. Conversely, models where, for example, a whole pentamer or hexamer is considered to be a lumped element, may not be viable either (i.e., the internal dynamics of these entities may be a key element of self-assembly as they deform during assembly). Furthermore, these internal degrees of freedom provide a sink for friction-associated energy transfer. The Smoluchowski equations we obtain are correct to even though the Liouville equation need only be solved to . Thus, our procedure is ideal for deriving higher order (augmented) Smoluchowski equations even for complex systems involving multiple OPs. This feature of our workflow follows from the use of the general equation for the reduced probability density that follows directly from the Liouville equation. These results enable the algorithm for BNS simulations suggested in
Fig. 2 . The Smoluchowski equation is solved in a Monte Carlo fashion via the equivalent Langevin equation with diffusion coefficients and averaged forces evolving with the OPs. While the MD timestep must be shorter than the 10 -14 second timescale of fast atomic vibrations and collisions, those for the present algorithm are limited by the characteristic timescale of the phenomena of interest (e.g., 10 -3 seconds for viral structural transitions). The resulting algorithm can therefore be many orders of magnitude faster than MD, despite the overhead from the computation of average forces and diffusion coefficients. The computation of diffusion coefficients must be carried out by constructing a time correlation function with evolution generated by the operator , and not by Newton’s equations directly. Thus, one may not use available MD codes to construct them. However, as the OPs are slowly varying, standard MD codes can be used to compute the diffusion coefficients approximately as long as the MD code is only run for a time short compared to the characteristic time of OP evolution, i.e., the constancy of the OPs must be checked during the construction of the diffusions coefficients. 17 Acknowledgements
This project at the Center of Cell and Virus Theory was supported by the U.S. Air Force, the U.S. Department of Energy, the Lilly Endowment, Inc., the Oak Ridge Institute for Science and Education, and the National Institute of Health.
Literature Cited
1. Shelley JC, Shelley MY, Reeder RC, Bandyopadhyay S, Klein ML. (2001).
Journal of Physical Chemistry B . : 4464 2. Murtola T, Falck E, Patra M, Karttunen M, Vattulainena I. (2004). Journal of Chemical Physics . : 9156 3. Izvekov S, Voth GA. (2005). A Multiscale Coarse-Graining Method for Biomolecular Systems. Journal of Physical Chemistry B . (7): 2469-2473 4. Izvekov S, Voth GA. (2005). Multiscale coarse graining of liquid-state systems. Journal of Chemical Physics . : 134105 5. Liu P, Izvekov S, Voth GA. (2007). Multiscale Coarse-Graining of Monosaccharides. Journal of Physical Chemistry B . (39): 11566-11575 6. Zhou J, Thorpe IF, Izvekov S, Voth GA. Coarse-Grained Peptide Modeling Using a Systematic Multiscale Approach. Biophysical Journal . (12): 4289-4303 7. Speelman, B., B.R. Brooks and C.B. Post (2001). Molecular Dynamics Simulations of Human Rhinovirus and an Antiviral Compound. Biophysical Journal , 80, 121-129. 8. Ortoleva, P. (2005). Nanoparticle Dynamics: A Multiscale Analysis of the Liouville Equation.
J. Phys. Chem ., 109, 21258-21266. 9. Miao, Y. and P. Ortoleva (2006a). All-atom multiscaling and new ensembles for dynamical nanoparticles.
J. Chem. Phys.
J. Chem. Phys.
Ap. J.
J. Chem. Phys.
Physics Letters
A69(5):367-369. 14. Bose, S., S. Bose and P. Ortoleva (1980). Dynamic Padé approximants for chemical center waves.
J. Chem. Phys ., , 4258-4263. 15. Bose, S., M. Medina-Noyola and P. Ortoleva (1981). Third body effects on reactions in liquids. J. Chem. Phys.
Faraday Discuss. Chem. Soc . , 1-20. 17. Ortoleva, P. (1992). Nonlinear Chemical Waves; John Wiley and Sons: New York. 18. Shea, J. E. and I. Oppenheim (1996). Fokker-Planck equation and Langevin equation for one Brownian particle in a nonequilibrium bath. J. Phys. Chem.
100 (49): 19035-19042. 19. Shea, J. E. and I. Oppenheim (1997). Fokker-Planck equation and non-linear hydrodynamic equations of a system of several Brownian particles in a non-equilibrium bath.
Pysical A
247 (1-4): 417-443. 20. Peters, M.H. (1998). Fokker-Planck equation and the grand molecular friction tensor for combined translational and rotational motions of structured Brownian particles near structures surface.
J. Chem. Phys ., Vol. 110, No. 1, pp. 528-538. 18 21. Peters, M.H. (1999). Fokker-Planck equation, molecular friction, and molecular dynamics for Brownian particle transport near external solid surfaces.
J. Statistical Phys ., Vol. 94, Nos. 3-4, pp. 557-586. 22. Coffey, W. T., Y.P. Kalmykov and J.T. Waldron (2004). The Langevin Equation With Applications to Stochastic Problems in Physics, Chemistry and Electrical Engineering.
World Scientific Publishing Co .: River Edge, NJ. 23. Jaynes, E.T. (1957).
Physical Review,
J. Phys. Chem. , Physical Review A . 136 (2A): A405-A411 30. Alder BJ, Wainwright TE. (1967). Velocity Autocorrelations for Hard Spheres.
Physical Review Letters .; 18 (23): 988-990 31. Alder BJ, Wainwright TE. (1970). Decay of Velocity Autocorrelation Function.
Physical Review A . 1 (1): 18-21 32. Ernst MH, Hauge EH, van Leeuwen JM. (1970). Asymptotic Time Behavior of Correlation Functions.
Physical Review Letters . 25 (18): 1254-1256 33. Zwanzig R, Bixon M. (1970). Hydrodynamic Theory of the Velocity Correlation Function.
Physical Review A . 2 (5): 2005-2012 34. Dorfman JR, Cohen EGD. (1970). Velocity Correlation Functions in Two and Three Dimensions.
Physical Review Letters . 25 (18): 1257-1260 35. Boon JP, Bouiller A. (1976). Experimental observations of long-time-tails.
Physics Letters A . 55 (7): 391-392 36. Paul GL, Pusey PN. (1981)Observation of a long-time tail in Brownian motion.
Journal of Physics A – Mathematical and General . 14 (12): 3301-3327 37. Weitz DA, Pine DJ, Pusey PN, Tough RJA. Nondiffusive Brownian Motion Studied by Diffusing-Wave Spectroscopy.
Physical Review Letters . 1989; 63 (16): 1747-1750 38. Zhu JX, Durian DT, Müller T, Weitz DA, Pine DJ. (1992). Scaling of transient hydrodynamic interactions in concentrated suspensions.
Physical Review Letters . 68 (16): 2559-2562 39. Kao MH, Yodh AG, Pine DJ. (1993). Observation of Brownian-Motion on the Time Scale of Hydrodynamic Interactions.
Physical Review Letters .70 (2): 242-245 40. Morkel C, Gronemeyer C, Glaser W, Bosse J. (1987)Experimental Evidence for the Long-time Decay of the Velocity Autocorrelation in Liquid Sodium.
Physical Review Letters .; 58 (18): 1873-1876 41. Vladimirsky V, Terletzky Y. (1945). Hydrodynamical theory of translational Brownian motion.
Zhurnal Eksperimentalnoi I Theoreticheskoi Fiziki .; 15: 258-263 /in Russian/ 19 42. Lisy V, Tothova J. (2004). On the (hydrodynamic) memory in the theory of Brownian motion. arXiv:cond-mat/0410222v1 . 43. Jaqaman, K. and P. Ortoleva. (2002). New space warping method for the simulation of large- scale macromolecular conformational changes.
J. Comp. Chem.
23: 484-491. 44. Miao and Ortoleva (2008). Molecular Dynamics/OP eXtrapolation (MD/OPX) for Bionanosystem Simulations. J. Comp. Chem. (submitted).
45. Shreif, Z. and P. Ortoleva. (2007). Curvilinear All-Atom Multiscale (CAM) Theory of Macromolecular Dynamics.
J. Stat. Phys. (in press). 46. Shreif and Ortoleva (2008). Computer-aided design of nanocapsules for therapeutic delivery Computational and Mathematical Methods in Medicine. (submitted). 47. Shreif and Ortoleva (2008). Derivation of Augmented Smoluchowski Equations Through Higher Order Multiscale Perturbation Analysis. (in preparation). 48. Pankavich, S., Z. Shreif, and P. Ortoleva (2008). Multiscaling for Classical Nanosystems : Derivation of Smoluchowski and Fokker-Planck Equations.
Physica A . (to appear) 20
Appendix A: State Counting
The first step in ensemble construction as in Section II is to develop a way to estimate the number of distinct quantum states in the 6 N dimensional volume element of atomic position/momentum space. For example, there is one state in a volume for N identical atoms. If is a product of M +1 Dirac delta functions ( M for the structural OPs and one for the system energy), then the aforementioned state count is given by where specifies the state of the N -atom system over which integration is taken. The state density factor is constructed as follows. Let be the unit step function: . Then is 1 for and zero otherwise. If is small, then . But for Dirac delta function . With this, . The factor is zero except in a small zone of values around of volume for width in each of the OPs, . Let be the product of the state counting factors for each of the types of atoms in the system. Then, the number of quantum states in the zone for which the product of D factors is one, is given by . Comparing this result with ( A.1 ) for small implies . The are all small constant quantities. Thus, is a constant, i.e., independent of . With this counting of states for the iso-OP systems, in Section II we make information theoretic arguments to construct for a range of distinct experimental conditions relevant for the study of BNS as follows.
Appendix B: Fourier Transform Method
In our earlier study, S was first maximized with respect to , constrained by normalization and the ensemble average values of a set of OPs . In the present notation, this becomes fixed. ( B.1 ) Introduce Lagrange multipliers associated with the M +1 constraints ( B.1 ) for the M +1 OPs constrained entropy maximization. One obtains the conditional probability given by ( B.2 ) . number of quantum states available to the system with in a fixed, narrow range about , ( A.1 ) 21 For the second level, we develop a more general solution in the form of a linear combination of the with unit-normalized weight ; in analogy with our earlier study , we find that the more general distribution, denoted , takes the form ( B.3 ) (
B.4 ) for ( M +1)-fold delta function and partition function and coarse-grained probability . The distribution ( B.3 ) is microcanonical in character for the detailed atomic configurations consisting of all states with fixed value of the OPs; thus, it has an equi-a priori principle character, i.e., all states with a given are equally likely in the absence of additional information. The result (
B.3) is derived from (
B.2 ) using properties of the Fourier transform. Consider a function of a set of L variables. The transform for the set of L “wave vectors” k is defined via ( B.5 ) The inverse relation reads , (
B.6 ) where . For functions and one has . ( B.7 ) The result (
B.7 ) follows from the inversion formula which yields ( B.8 ) Thus, is related to the integration of / N over all imaginary OP values. As is arbitrary, further assumptions related to the nature of that stem from the problem of interest must be introduced, e.g., the initial data. Finally, the thermal case with corresponds to keeping the energy fixed only by its average so that . Appendix C: Constructing order parameters from orthonormal polynomials
Consider a nanostructure embedded in a volume . Basis functions for a triplet of labeling indices are introduced which are orthonormalized. In this method, the nanosystem 22 deforms in 3-D space, which is considered to be a displacement of an original point . Deformation of space taking any to is continuous and is used to introduce OPs via . (C.1)
As the change, space is deformed, and so does the nanosystem embedded in it. The objective is to ensure that the dynamics of the reflects the physics of the BNS and that the deformation reflects key aspects of the atomic-scale details of the structure. In this way, the constitute a set of vector OPs that serve as the starting point of our multiscale approach. In our approach, atom is moved from its original position via the above deformation by evolving the and correcting for atomic-scale details as follows. Given a finite truncation of the sum in (C.1) , there will be some residual displacement for individual atoms. Denoting this residual for atom as . (C.2)
The size of can be minimized by the choice of basis functions and the number of terms in the sum. Conversely, imposing a permissible size threshold for the residuals allows one to determine the number of terms to include in the sum. To start the multiscale analysis, the must be expressed in terms of the fundamental variables . Let be the mass of atom . Multiplying (C.2) by and summing over the N atoms in the system, one obtains ( C.3)
The integral orthogonality of the basis functions implies that the matrix is nearly diagonal. Thus, the OPs can easily be computed in terms of the atomic positions by solving (
C.3) numerically. More specifically, when most of the system is occupied with atoms, the i sum is essentially a Monte Carlo integration. The orthonormality of the basis functions implies that and ( C.3 ) can be approximated as (
C.4 ) The contribution is neglected in arriving at this definition of as fluctuates in direction and magnitude with , while the basis functions that capture overall characteristics of the nanosystem (e.g., position, orientation, nanoscale structure, and hydrodynamic modes) vary smoothly by design. Thus, (C.3) is not an approximation; rather, the above is a way to argue for the definition of OPs that express coherent behaviors of the nanosystem. With this definition of the , (C.3) is an exact relationship, since the correct errors in the displaced atomic positions over-and-above the coherent contribution from the sum.
Our approach has several conceptual and technical advantages. In the above we include all atoms in the system, e.g., those in the nanostructure and its microenvironment. This captures the boundary layer of water that tends to accompany a nanoparticle due to the viscous nature of water at the nanoscale and the interaction forces. Thus, the method captures hydrodynamic modes in the host fluid and layers of water bound to the nanoparticles and other structures. Inclusion of in the above expressions gives the character of generalized, center-of-mass (CM) variables. For example, if is a constant for a specific , the corresponding value of is proportional to the CM of the system. As for momentum of atom i , one has ( C.5 ) While has a sum of N atoms, many of which have similar directions due to the smooth variation of with respect to , the momenta have fluctuating direction and tend to cancel near equilibrium. Hence the thermal average of is small, tend to evolve slowly, and the ratio of the characteristic time of to that of atomic vibrations and collisions should be on the order of the number of atoms in the system, i.e., O( N ). This suggests that the are slowly varying, and therefore satisfies a key criterion to be an order parameter and serve as the starting point of our multiscale analysis. As N is large, it is convenient to define the smallness parameter (where is the mass of the nanosystem and mm