Methods and Models

The most direct method to simulate protocellular functions at a molecular level is molecular dynamics [13]. In this method, Newton's equations of motion for all the atoms in the system are solved as a function of time. Exactly the same approach is commonly used to study computationally large systems in chemistry and structural biology.

The first task is to define a protobiologically relevant, yet computationally tractable, model system. As a suitable model for membrane forming material we selected glycerol monooleate (GMO). A GMO molecule is composed of a glycerol head group linked by an ester bond to a hydrocarbon tail containing 18 carbon atoms with a double bond in the middle. Undoubtedly, the GMO bilayer does not faithfully represent the composition of protocellular membranes which, most likely, were built of highly heterogeneous amphiphilic material that contained both charged and uncharged molecules. However, in contrast to GMO membranes, such heterogeneous systems cannot be reliably modeled by computational methods. Even though the GMO bilayer has a homogeneous composition, it still retains important features expected of primitive membranes, namely structurally simple head groups and a highly fluid interior. In these respects it may be a better protocellular model than membranes built of phospholipids that form the walls of contemporary cells. This point is further reinforced by difficulties in identifying sufficient sources of phosphate on the early Earth [14] to allow for the formation of significant quantities of phospholipids.

Present computational resources do not permit simulating a whole vesicle in an aqueous environment. Instead, we considered only a part of this system -- a fragment of the membrane consisting of 72 GMO molecules and covering an area of 37 37 Ų, surrounded by 2300 water molecules evenly distributed on both sides of the bilayer. Undesirable discontinuities at the edges of the system were eliminated by applying periodic boundary conditions. This system is sufficiently large to study such phenomena as transport, organization and chemical reactions of molecules at water-membrane interfaces, and the formation of transmembrane proton gradients.

The equations of motion describing the system were solved numerically on a step-by-step basis using a finite difference method. From dynamical information about the system at time t, we obtained the positions and velocities of all the atoms at time that, in turn, were used to calculate these quantities at , etc. This procedure was repeated many times, resulting in a complete, microscopic description of the system as a function of time (called a trajectory). Macroscopic quantities can also be obtained from this trajectory. In our simulations, was set between 210^-15 and 510^-15 sec and the length of trajectories varied between 10^-9 and 10^-8 sec. Thus, each trajectory required approximately 10⁶ time steps.

In each step, the forces acting on each atom in the system have to be calculated. These forces are derivatives of the potential energy function with respect to the atomic coordinates. The potential energy function was represented as a sum of contributions from electrostatic and van der Waals interatomic interactions as well as terms describing intramolecular bond and angle vibrations and changes in the dihedral angles formed by three consecutive bonds. Electrostatic contributions were evaluated as a sum of Coulomb energies between partial point charges assigned to atoms. For water, the TIP4P potential energy function [15] was used. Potential energy functions for GMO were developed by Wilson and Pohorille [16] and for peptides by Cornell et al. [17].

Despite large computational effort, the time scale covered by molecular dynamics simulations remains quite short. The probability of observing processes that typically occur at considerably longer time scales in such simulations is very low. However, reliable structural and energetic information about these processes can still be obtained by dividing them into several consecutive stages that are simulated separately. For example, solute transport across the membrane can be represented as a series of stages in which the solute is progressively moved across the membrane in the direction perpendicular to the bilayer (the z-direction). At each stage, the solute is constrained to lie within a ``window'' along z. The free energy change at each stage, , as a function of the parameter z can be obtained by observing the probability, P(z), of finding the solute at z:

 

where k_b is the Boltzmann constant and T is the temperature of the system. If the ranges of z explored by the solute in consecutive stages overlap, the dependence of on z for the whole process can be obtained from the condition that the free energy must be a continuous function of z. Exactly the same procedure can be used to study other processes of interest described by different parameters.