1. Introduction

An ongoing challenge in the field of biological computer simulation involves the treatment of disparate time and length scales. Many bioassemblies (for example, the DNA double helix) explicitly contain both microscopic and macro-scale dimensions in their structure. The length along the strand of the molecule extends to macroscopic dimensions, while its width is microscopic. Thus, one can think of the “width” of the DNA strand as existing in a microscopic temporal/spatial regime, whereas the “length” is macro-scale in nature. Other biostructures such as cell membranes exhibit similar characteristics.

The characteristic of multiple, disparate time and length scales in the biological structures of living systems presents a clear barrier to their computational modeling. To properly model trans-temporal and trans-spatial structures involves crossing this microscopic/macroscopic barrier to incorporate an adequate degree of information transfer between these disparate scales.¹

On one hand, simulations on the scale of millions of atoms, taking months of computational time, do not begin to approach the time scales required to examine the temporal interplay between atomic- and cellular-level processes. Even advanced simulation techniques employing multiple time-step dynamics are still, at a fundamental level, attempting to describe macroscopic-level phenomena with microscopic-scale models. Conversely, the most detailed nonlinear equation cannot explicitly model molecular-level interactions if the information from such interactions is not properly incorporated.

In the context of bio-assemblies and computer simulation, the capability to interlink atomistic detail with macroscopic events must involve a fundamental new direction in simulation methodology. An intrinsic coupling between these disparate length and time scales that retains inherent nonlinear behavior but avoids the computational trap of unrealistically large simulations is required.

Molecular dynamics simulations can model systems of the order of 10⁵ atoms for time scales of the order of nanoseconds, which is still orders of magnitude below continuum length scales and time scales, these in turn being of the order of millimeters and milliseconds. Continuum-level simulations operate in the latter regimes, but do so at the expense of explicit molecular-level detail and structure. Instead, the physical properties of the material are contained within various transport coefficients which appear in the constitutive relations that are a crucial step in the formulation of the continuum-level dynamics. Generally, these transport coefficients (which could include shear or bulk viscosity, Young's modulus, or Poisson's ratio, for example) are considered as initial input to the simulation.

However, at some fundamental level, all transport coefficients and equilibrium properties such as pressure and internal energy are averages of complex microscopic-scale interactions. Thus, one might already have deduced that there is an obvious “link” between micro- and macro-scale dynamics that is inherently embedded within the transport coefficients themselves. If the relevant transport coefficients required for a continuum-level simulation are calculated from a detailed atomistic model, we can essentially “jump” from the atomic spatial/temporal regime to the macro scale.

The converse is also true. The parameters that define a microscopic system are quantities such as density, temperature, and pressure. These parameters, like the transport coefficients in the continuum case, are usually considered as initial input. Moreover, quantities such as the density are inherently macroscopic in nature. Thus, if the density for an atomistic-scale simulation is obtained from the calculated density in a corresponding continuum-level simulation, we have essentially performed the reverse time/space “jump.”

The combination of these two input/output sets can result in a closed feedback loop between micro-scale and macro-scale temporal and spatial regimes. In many cases the approximation that transport coefficients (inherently microscopically determined quantities) are not altered by continuum-level fluctuations in density is reasonable. However, there are a number of systems in which this approximation is not valid, and some form of feedback or information transfer across these disparate time scales is required to completely elucidate the dynamics of the system. For example, highly non-elastic materials could have density-dependent transport coefficients; that is, the value of the required transport coefficient in the constitutive relation will depend on the density of the material at that region. In a situation such as this, a method of calculating the density-dependent transport coefficient at that region is required.

For this paper, we present a new interface between microscopic- and macroscopic-level simulations. The general idea, as previously mentioned, is to calculate the transport coefficients required for a continuum-level simulation from corresponding molecular dynamics simulations. However, the interlink between the two simulations allows the examination of the dynamic feedback between microscopic and macro-scale time and length regimes. The material to be studied is represented both at the continuum and microscopic levels, and selected information is transferred between the two levels. The present study is restricted to biological membranes, as these systems represent an interesting template for the examination of trans-temporal and trans-spatial phenomena. In their inherent nature, one direction (the membrane width) exists in a microscopic domain, while in perpendicular dimensions the length scales are effectively macroscopic relative to the membrane thickness. It should be noted that the technique, which is referred to as “dynamic feedback,” is completely general and can be applied to a wide array of systems.

The paper is divided into several sections. In Section 2 we discuss some general properties of real biological membranes, and how they are commonly modeled with computer simulation. We then present in Section 3 a brief overview of the microscopic simulation technology that is employed; non-equilibrium molecular dynamics (NEMD) is used extensively as a means of calculating transport coefficients from microscopic models. Aspects and details of the continuum-based simulation technique, known as the material point method, are discussed next, in Section 4. The dynamic feedback mechanism is described in Section 5. Finally, we present results in Section 6 for the dynamic feedback between the microscopic and continuum levels for a spring-membrane experiment, along with the future outlook for the method in Section 7. Concluding remarks are given in the last section.

2. Biological membranes

The cell membrane of all living organisms defines the boundary between the outside environment and the biological inner system. Furthermore, it not only defines the boundary of the system, but also regulates the flux of materials entering and exiting the cell. The thickness of the lipid membrane is microscopic; however, its surface area extends to macroscopic dimensions relative to its thickness. At an atomic level, the membrane is composed of lipid molecules, which contain a hydrophilic (water-attracting) group at one end, and then two hydrophobic (water-repelling or “oily”) carbon chains. Thus, the structure of the lipid molecule could be described generally as elongated, much like an ellipsoid of revolution or a spherocylinder, in which one end is attracted to water and the rest is repelled by it. At specific concentrations and temperatures, a lipid bilayer is formed. In this phase, the lipid molecules arrange themselves side by side such that their hydrophilic ends point outward to the surrounding solvent. The membrane is composed of two lipid layers, with the hydrophobic (oily) components pointing inward and the hydrophilic ends in contact with the solvent. Moreover, the lipid molecules themselves diffuse within the plane defined by the membrane bilayers; thus, the membrane is best described as a liquid crystal, or gel, rather than a solid. The fluid-like nature of the lipids within the plane of the membrane is absolutely critical for cellular functions such as ion transport.

There are various models for biological membranes ranging from very detailed atomistic models of lipid bilayer systems [1, 2] to greatly simplified models (hard spherocylinders, ellipsoids of revolution) used in liquid crystal studies [3, 4]. For this study we chose an amorphous membrane model that can easily be constructed and simulated. A snapshot of the microscopic model of the membrane is shown in Figure 1. Basically, the membrane is constructed from a random configuration of spherical atoms that are bonded together by harmonic spring bonds. The bonds between the atoms ensure that the atoms remain bound within the membrane, and the resulting material can be thought of as a toy “cross-linked polymer.”

Figure 1

3. Non-equilibrium molecular dynamics

Classical equilibrium molecular dynamics (MD) typically involves solving Newton's equations of motion for a relatively small number (usually N < 10 000) of atoms and molecules. If Newton's equations of motion for the particles are solved, F = ma, where F is the force acting on a particle, m is the mass, and a is the acceleration, the total energy E of the system is conserved. However, quite often one wishes to constrain the temperature T or the pressure P to be constant. It is possible to generate “synthetic” dynamics [5] in which the equations of motion are extended such that the new dynamics conserve T or P rather than E. In the large-system limit, where the fluctuations in these quantities become infinitesimally small, these synthetic equations of motion approach Newton's equations.

Molecular dynamics can also be further extended to the non-equilibrium regime. With non-equilibrium molecular dynamics (NEMD), an external field, force, or flux is imposed on the system [5–7]. The non-equilibrium state produces heat because of the work being performed on the system. The same synthetic dynamics used to maintain constant temperature or pressure can be employed to remove the excess heat. If the heat produced can be removed, the system can, under some circumstances, approach a non-equilibrium steady state. The subtle aspects of the statistical mechanics of NEMD are discussed in detail in Reference [5], and they are not discussed here. For our purposes, NEMD is an efficient and elegant method to calculate transport coefficients, in which the transport coefficient of interest is calculated in the limit that the non-equilibrium perturbation goes to zero. The advantage of this approach over equilibrium techniques for the calculation of transport coefficients is that NEMD does not usually require either long or large simulations. The efficiency of this approach for atomistic membrane simulations is a key feature of this paper and is an important part of the bridging strategy between the micro- and macro-scales.

Cyclic compression NEMD
Cyclic compression NEMD has been used to calculate the bulk viscosity of fluids, and has been very successful in comparison with traditional methods (Chapman–Enskog) [8, 9]. This method can be extended to viscoelastic solids, as well as to materials with non-elastic properties. The idea is to introduce an artificial volume oscillation and then to correspondingly calculate the force required to generate the oscillation. For example, if an elastic material is compressed in some direction, the force required for that compression occurs in the same direction, and the material responds with a force in the opposite direction. When the compression is in the x direction, a relationship between the oscillation and the force required for its generation (in the case of an elastic material) is

x=E,

(1)

where _x is the force per unit area in the x direction, = x/L is the dimensionless dislocation, and E is the modulus. If E is time-independent, the time derivative of this equation relates how the force changes are related to the oscillation changes,

x =E ,

(2)

where = u_x/x and u_x is the flow field in the x direction. In NEMD, it is convenient to use the pressure tensor P rather than the external force per unit area , where P = –; thus, the relationship between the time derivative of the pressure component and the strain rate is

–xx=E.

(3)

The modulus is calculated in the limit that the frequency and amplitude of the oscillation approach zero. With cyclic compression NEMD, we impose a dilation rate given by
= cos(t), where is the amplitude of the oscillation (dimensionless), = 2/, and is the wavelength.

4. Continuum-level simulation and the material point method

The material point method (MPM) is a numerical technique for solving large-deformation problems in continuum mechanics [10–12]. In MPM, a fluid or solid body is discretized using an unconnected set of material points that are followed throughout the deformation history. Details of the continuum-level equations can be found in Reference [10]; here we describe only the technique. The idea with MPM is that a material is modeled at the macro-scale level by treating it as a smooth, continuous system. There is no molecular-level granularity whatsoever; however, for computational purposes the material is divided into a number of material “points” that have mass, velocity, location, and stress. The material points move according to the forces acting on them, with the forces being determined by the collective movements of the other material points, as well as any external stresses.

The basic equation of motion solved via MPM is given by

• + b = a,

(4)

where is the stress tensor, b is the specific body force, is the mass density, and a = is the acceleration. The actual algorithm involves “mapping” quantities such as the material point mass, the velocity, and the acceleration from the material points themselves onto a lab-fixed set of grid nodes. The grid nodes do not move. This transformation from the material point representation to the grid node representation is necessary in order to calculate the forces acting on the material points. Because the MPM method is completely general, it can essentially model any material at the continuum level. In the case of a biological membrane as previously described, solution of these dynamics requires a constitutive relation such as those given in Equation (2) and Equation (3). Once the modulus is known, the value of _x or –P_xx can be found from , and then, in turn, the acceleration can be calculated. The mapping of properties such as mass and momenta from the material points to the lab-fixed grid enables complete spatial freedom of material points. For dense fluids, liquid crystals, or gels, this extension is necessary, since a material point can essentially move about the entire region.

5. Dynamic feedback between MD and MPM

The NEMD/MPM interface can be represented diagrammatically as shown in Figure 2. Beginning with the topmost image (1), the NEMD/MPM interface is initialized with NEMD simulations of the membrane at a microscopic level, where initial transport coefficients (elastic moduli or viscosities) are calculated. With these initial values, an MPM simulation of the membrane is begun (2). As the MPM simulation evolves in macroscopic time, new densities are computed and used as input in a subsequent set of microscopic simulations. On the microscopic scale, the system has now visited higher or lower densities (3), and new NEMD calculations at these new states are performed. New transport coefficients are then fed back into a new MPM simulation (4). This iterative process continues in a loop, and the inherent nonlinear and complex behavior is ultimately incorporated. Small instabilities at either the micro or macro level can cross-propagate temporal regimes, resulting in membrane deformation that could lead, for example, to membrane rupture. At a deeper level, what in essence results from this micro-to-macro feedback is that Equations (2) and (3) explicitly contain E[(r(t), t)], which is a nonlinear equation in space and time (to be contrasted, for example, to a perfect elastic solid).

Figure 2

At the microscopic level, as long as the temperature and potential function are not altered, the modulus will depend only on the density, i.e., E[(r(t), t)]. Since the value of E is unique to the density (as long as the model and temperature at the microscopic level are not altered), rescanning previously sampled densities with NEMD is redundant. In this case, the iterative feedback loop as previously described can be significantly simplified by constructing a table of E versus . This table can also be constructed “on the fly”: After different densities have been sampled by the MPM simulation, a series of NEMD simulations are initiated, and new values of E corresponding to the sampled densities are found. For the system studied in this paper, the dynamic feedback mechanism takes the form of a precalculated table of E versus , which is sufficient as long as the preliminary set of NEMD simulations to determine the density dependence of E spans the accessible densities sampled by MPM.

6. Results

NEMD simulation
The model membrane for this study is as described in Section 2. For this system, we employed “molecules” of diameter = 5 Å, mass m = 15 amu, at a temperature of T = 298 K, which reasonably models an atomistic-level membrane. Figure 1 is a snapshot of the system, in which the random bonds between the atoms are clearly shown. Briefly, the NEMD simulations were performed by first running an equilibrium simulation of N = 900 molecules for a time of 100 ps, and then the system was subjected to a series of 40 to 100 non-equilibrium cyclic volume compressions in the x-direction of amplitude and frequency .

The density dependence of the modulus E was found by performing the series of NEMD cyclic compression runs at five different mass densities: = mN/V = 0.096 amuÅ^–3, 0.087 amuÅ^–3, 0.080 amuÅ^–3, 0.049 amuÅ^–3, and 0.024 amuÅ^–3. (Note that = 0.12 amuÅ^–3 corresponds to a simple cubic close-packed structure.) For each density, two different cyclic compression amplitudes corresponding to = 0.265 and 0.132, and five different frequencies ranging from = 0.07 ps^–1 up to = 0.29 ps^–1, were used. The magnitudes of the amplitudes were less than one percent of the original cell size, and the shortest frequency corresponded to a wavelength of 50000t, where t = 0.001 ps is the time step used in the simulation. For this system, this wavelength is almost an order of magnitude longer than the typical time necessary for the system to relax to equilibrium from an initial configuration. In Figure 3 the modulus is plotted versus dilation frequency for five different densities. Except for the most dense fluid system at = 0.096 amuÅ^–3, E is basically independent of the compression frequency, and the extrapolated value is not significantly different from the non-equilibrium values. At = 0.096 amuÅ^–3, a small increase in E is observed for increased frequency. Results for two different amplitudes are shown (solid and open squares), and the extrapolated value of E is essentially independent of amplitude.

Figure 3

As an alternative method of calculating E, and to check whether or not there could be very low-frequency or low-amplitude nonlinear behavior, an elastic expansion NEMD simulation series at the highest density ( = mN/V = 0.096 amuÅ^–3) was performed. Very briefly, the elastic expansion technique rigorously should be applied only to elastic materials in which the stress versus strain curve is linear. An external stress is imposed on the system, and the system is allowed to expand in one direction, until a new “steady state” dislocation is reached. This amounts to stretching the material and calculating the dimensionless displacement that results from the stretch. With this technique we find excellent agreement with the calculated value of the modulus at = 0.096 amuÅ^–3, with E = 1.126 × 10⁵ kgm^–1s^–2.

The extrapolated values of E can be plotted versus = Nm/V. In Figure 4 the nonlinearity of E with respect to density is clearly shown. At low densities, E remains almost constant with (the “harmonic regime”), while at fluid densities with > 0.07 amuÅ^–3, highly nonlinear behavior is observed. The nonlinearity of E at fluid densities arises from atoms being pushed together with strong short-range repulsion interactions; thus, the material rapidly becomes “stiffer” at densities beyond > 0.07 amuÅ^–3.

Figure 4

MPM simulation
The results from Figure 4 (in which a fourth-order polynomial was used to fit the data) were then employed as initial input in an MPM simulation of the membrane. Recall that in a traditional MPM simulation, the modulus would be considered as initial input to the program, independent of the density of the material. In this simulation, the modulus as a function of the density is the initial input, which has been uniquely determined for the chosen membrane model, at a specified temperature.

The MPM simulation was set up to model a membrane-spring mass experiment in which a macroscopic membrane is anchored at one end and a weight attached to the other. Under the force of gravity, the weight bounces. For this experiment, a membrane of length 0.3 mm and cross-sectional area 0.001 mm² was chosen. Ten material points and four grid nodes were placed on the membrane, and the tenth material point (the one at the bottom) was given a mass 10000 times greater than the other material points to model the weight attached to the end of the membrane. The initial value of the modulus was 1.126 × 10⁵ kgm^–1s^–2, corresponding to a mass density of = 0.096 amuÅ^–3. Further details are summarized in Table 1, and a diagram of the continuum-level model is shown in Figure 5. This particular model is presented here as a demonstration of the NEMD/MPM feedback loop. Other, more realistic, simulations can be constructed, but for the present purposes this continuum problem is attractive in that it has an analytic solution in the density-independent case.

Figure 5

Figure 6 presents the results for the dynamic feedback interface of a microscopic-level simulation with the continuum-level MPM model. The curves in the plot track the displacement of the “heavy” material point over time. Note that the spatial displacement is now given in millimeters, while the time is in milliseconds. Recall that the modulus was calculated using a microscopic-scale simulation on a length scale of the order of nanometers and time scales of the order of picoseconds. An informative calculation is to determine the system size and the computational time required to model this simple spring experiment using a detailed microscopic model with no bridging to the continuum-level MPM simulation. To match the chosen dimensions of the MPM simulation, a microscopic simulation of 1.9 × 10¹⁵ molecules running for 1.0 × 10¹⁰ ps would be required. Present atomic-level simulation capabilities are nowhere near being able to handle such simulations. In contrast, with NEMD/MPM feedback, the MPM calculation takes about a minute on a workstation (SGI 02, R5000, 200 MHz), while each NEMD modulus point takes about five minutes to calculate.

Figure 6

In the present application, the actual implementation of the dynamic feedback is contained within the MPM algorithm, in which the density about each material point is integrated in time. When the density at a specific material point is updated in the time integration, the value of E, as shown in Figure 4, is employed in the constitutive relation

xx(t)=–E[(t)]t).

(5)

As the membrane extends, the density within it decreases, and from Figure 4 it can be seen that E decreases with in the sampled density regime. Figure 6 shows that the membrane extends slightly farther when dynamic feedback is included (solid curve) than in the constant-E case (dotted curve).²

7. Future applications to biological membranes

We have demonstrated that it is possible to interface microscopic- and macroscopic-level simulations by constructing a dynamic feedback loop within an NEMD/MPM simulation. The key to the method lies in identifying those microscopic quantities that must be transferred to the macro scale and, also, which macroscopic properties are propagated backward to the microscopic regime. As a first demonstration of this technique, we have chosen a simplified model that not only demonstrates the basic algorithm but also clearly elucidates the effects of including dynamic feedback.

For the future, there are a number of problems in computational biology where our methodology will help overcome the computational barrier to long time- and length-scale simulations. In the context of membrane simulations, extending our approach to study these systems will involve basically two steps:

• Employing a more detailed microscopic model of the lipid membrane as previously described in Section 2.
• Extending the MPM model to describe membranes in more realistic geometries. For example, the MPM simulation would model a membrane in a “drum” geometry, or even a “bubble” geometry.

There are currently two problems of immediate interest that will inherently require some sort of micro-to-macro feedback to completely model the full scope of their behavior:

1. Osmotic effects in black membranes
Experimentally, synthetic lipid bilayers, or black membranes, are used as a model system for real cell membranes. Black membranes can be formed by creating a small hole (~0.8 mm) on a Teflon barrier that is immersed in water. The membrane is formed by introducing the lipid molecules at the small hole between the two aqueous compartments on either side of the barrier. The inherent selectivity of biological membranes gives rise to osmotic pressure effects: If a solution of water and a substance impermeable to the membrane are placed in one compartment and pure water is in the other, a flux of water from the pure side to the solution side is observed. On a microscopic scale, water molecules percolate through the lipid membrane, driven by the chemical potential gradient. On a macroscopic scale, the membrane is subject to a pressure head directed from the solvent to the solution. The flux of water through the membrane results in a deformation of the membrane across the pore, and the deformation, manifested microscopically as a local change in lipid density, affects the flux. At some critical value of the osmotic pressure gradient, the flux will become too great, and the membrane will rupture.

To study this problem with NEMD/MPM feedback requires the extension of the MPM membrane simulation to model a two-dimensional membrane “disc” or “drum.” With the extension of MPM, a concurrent extension to the NEMD component is also required. For this model membrane two moduli are required: the shear and the bulk elastic moduli (or viscosities).

2. Nonlinear effects of biologically relevant biomolecules in cell membranes
The presence of nonlipid biomolecules such as cholesterol in cell membranes should alter the membrane's macroscopic material properties. This perturbed mechanical state should then feed back to the microscopic regime. There is already evidence from detailed microscopic-level simulations that the presence of nonlipid molecules can alter the local material properties of a membrane [1, 2]. That is, if a cholesterol “density” is present, the membrane generally becomes stiffer. However, from this result, without some sort of micro-to-macro interface, there is no way of resolving how these small changes in material properties will alter macro time- and length-scale dynamics.

The treatment of this problem will require incorporating more detailed and realistic potentials to model the lipid bilayer and the foreign “agent.” Systems of the scale of di-palmitoylphosphatidylcholine, or DPPC [1, 2], can be implemented, and the inclusion of more complex structures within the membrane, such as membrane-bound proteins, ion channels, and receptor sites [13, 14], can be incorporated to enrich the model within the NEMD/MPM framework. One possible scenario involves the inclusion of an ion channel within the lipid membrane. The NEMD/MPM feedback allows the inherent nonlinear aspects of the problem to be captured. If an unstable, or statistically improbable, state point is sampled at the continuum level, this information is transferred to the microscopic simulation. As the feedback continues, this nonlinear fluctuation has the ability to propagate in the macroscopic dynamics.

The present NEMD technique of cyclic compression can also be used in nonhomogeneous systems such as lipid bilayers. However, since this is a boundary-condition algorithm (that is, the non-equilibrium perturbation is introduced by altering the shape of the periodic lattice), all molecules within the system are inherently subject to the perturbation. Alternative methods analogous to the sinusoidal transverse force method [15] employ synthetic external fields, rather than altered boundary conditions, to drive the system from equilibrium. With these techniques, the external field acts on only those molecules of interest; thus, the mechanical properties of a lipid bilayer can be found even if the membrane is surrounded by a solvent. When NEMD is combined with hydrodynamic methods for evaluating the pressure tensor [16], its flexibility in highly inhomogeneous systems is greatly extended.

8. Conclusion

In this paper we have described a simulation methodology in which the computational barriers of disparate time and length scales that occur in many biological systems can be overcome. The method is general and relies on interlinking microscopic and continuum-level mechanics. In the context of biological structures, we have performed a series of computer experiments in which a membrane model is simulated in both the microscopic and continuum-level regimes. The necessary quantities required by the continuum-level simulations were calculated from the microscopic simulations, and in return the conditions for the microscopic simulations were defined by the continuum dynamics. The results from these computer experiments indicate that there may be subtle aspects of biological dynamics in which a dynamic information transfer between the atomic and the continuum level is necessary in order to fully model the behavior of the system.

Acknowledgment

This research was supported in part by the University of Utah Research Foundation Incentive Seed Grant Program.

References

Footnotes

Received June 20, 2000; accepted for publication August 8, 2000