format: text

text:
# Molecular dynamics simulation of polymer helix formation using rigid-link methods

## I Introduction

Polymers, because of their importance and complexity, have provided a longstanding challenge for computer simulation. Over the years, the field has become fragmented, both in terms of the problems addressed and the methodology employed. Broadly speaking, the kinds of system studied can be classified into distinct groups; there are biological heteropolymers, a category dominated by the proteins; homopolymers and block copolymers that include a great variety of molecular types, from alkanes to plastics; and idealized polymer models used for elucidating general principles such as the theta point, reptation, and multiphase behavior. The computational techniques span an equally broad range; they include molecular dynamics (MD) simulation employing models that represent the molecules at various levels of detail, ranging from fully atomic to highly reduced descriptions; Monte Carlo sampling of both continuum- and lattice-based systems, again with different levels of representation; and exact enumeration of small systems aimed at eliminating the sampling errors inherent in the other methods. While all three kinds of methodology provide important information about equilibrium behavior and, in a sense, amount to doing statistical mechanics numerically, it is only the MD approach that provides access to the dynamical and nonequilibrium aspects of the behavior; while it might be argued that Monte Carlo shares some of this capability, the associated dynamics is artificial and entirely a consequence of the chosen stochastic sampling algorithm, and so bears little relationship to Newtonian dynamics. Lattice-based approaches, though offering a vastly reduced configuration space, have the additional problem of the discreteness of the lattice on which the polymer is embedded, and the consequent absence of gradual transitions between different configurations.

The inherent difficulty in polymer simulation is that the problem naturally embraces a broad range of timescales, ranging from very fast processes associated with bond vibration, followed by the somewhat slower, highly localized conformational changes such as crankshaft motions, then the even slower aspects of reorganization such as the still relatively localized process of helix formation, and, finally, the typically extremely slow changes that lead to the emergence of tertiary structure characteristic of protein folding and to polymer diffusion in a concentrated solution. The timescales associated with this hierarchy of processes span a range considerably in excess of ten orders of magnitude, and so such systems are clearly not generally amenable to direct modeling, unless subjected to major simplification. Considerable effort has been invested in the design of models and simulation methods with the aim of alleviating this problem to at least some degree.

One especially important application of polymer simulation is in the field of protein folding, e.g., Brooks III et al. (1988); Dill et al. (1995); Shakhnovich (1997); Brooks III (1998); Pande et al. (2000); achieving an understanding of the mechanisms underlying this important process presents a major challenge to computational biochemistry. Protein modeling runs the gamut from, at one extreme, highly detailed molecular representations involving potentials derived from a mixture of theory and experiment, together with a solvent of individual water molecules, all solved by MD and an enormous amount of computational effort Duan and Kollman (1998); Berendsen (1998), through highly simplified models also solved by MD Levitt and Warshel (1975), to yet even simpler models embedded in lattices with only a limited number of degrees of freedom (DOFs) studied using a suitable Monte Carlo procedure and a greatly reduced investment in computing Dill et al. (1995); even complete enumeration of all conformations is sometimes feasible Crippen and Ohkubo (1998). While the manner in which the amino acid sequence of any given protein is able to determine its presumably unique spatial structure continues to be the subject of intense study, of no less importance is the question of the folding pathway – the preferred route (or routes) through multidimensional conformation space eventually terminating at the native state. While all the widely differing methodologies enumerated above can be used for studying folded states, the collective dynamical processes that underlie folding really demand an approach based on MD. But, after reaching this conclusion, there is a practical question of whether, even after substantial simplification, serious progress in understanding the mechanisms of folding can be achieved by computer simulation, owing to the diversity of intrinsic timescales; while substantial advances have been made, a great deal remains to be done before this question is answered.

The goal of the present paper is twofold. The first goal is a demonstration of a different perspective on the MD approach to studying protein folding. The most ambitious level of modeling is based on carefully constructed potential functions, often with a multitude of parameters; since the native conformation generally corresponds to the state of minimum free energy, establishing the details of these interatomic interactions, including solvent effects, provides the foundation for such work. Determining whether the known native state of a given protein is the one favored by energetic considerations is in itself a complex optimization task, but following the full dynamics over a sufficiently long period of time for the major structural changes that typify protein folding to occur verges on the impossible. The approach adopted here is just the opposite, and the question posed is the following: Given a known structural motif, such as the helix, and a simplified model of a polymer chain with a readily determined, unique ground state corresponding to this configuration, as in Fig. 1, will the chain collapse into this state within a reasonable amount of computation time when allowed to move freely in space, as shown in Fig. 2, while subjected to gradual cooling?

The most elementary of these organized structures is the helix, which, while being a prominent feature in many globular proteins, is only classified as a secondary structural element (the primary structure being the amino acid sequence itself), and because of its homogeneous nature (except for the ends) it might be argued that being able to fold a helix is not really a significant step in learning how to fold an entire protein. Therefore, another folding problem considered here is one with a ground state formed from an antiparallel pair of helices. This, too, is a recognizable element in some proteins, and is unquestionably classified as tertiary structure.

The obvious extension of this approach, a subject for future exploration, is to design simple models for other structural motifs, in the hope of learning more about folding by examining the collapse pathways of these idealized models; some structures might fold more readily than others, in which case the steric and topological issues involved could be investigated; for some structures there might be recognizable intermediate states along the folding trajectory; some cases might reveal useful properties that, when regarded as conformational (or reaction) “coordinates”, might serve in the design of other kinds of simplified models Crippen and Ohkubo (1998); and finally, once the simple version has been found to have the correct behavior, the models could be enhanced by gradually incorporating features from more realistic representations, including specific interactions and structural details. This represents the motivation for this kind of modeling approach.

The second goal is methodological. Even when considering the simplest of model polymers, in which, typically, all the molecular detail is absorbed into effective atoms located along the backbone chain (more so if this simplification is not made) there is a need to specify the internal DOFs of the system. One possibility is to assume that adjacent atoms are connected by stiff springs represented by a suitable potential function; in this case each atom has its full complement of three translational DOFs and, if these atoms are regarded as rigid particles rather than point masses, three rotational DOFs as well. If the bond potentials are made sufficiently stiff to correspond to a typical real system, the ensuing high-frequency vibrations impose a very small integration timestep, which runs contrary to the goal of efficiently simulating over long periods of time.

It is, however, possible to introduce geometrical restrictions, such as strictly constant bond lengths, while retaining a soluble dynamical problem. This is done by introducing holonomic constraints and Lagrange multipliers into the equations of motion Goldstein (1980), and then solving a set of algebraic equations while integrating the differential equations of motion. Two approaches have been developed for doing this; one involves initially solving the unconstrained equations of motion over a single timestep and then iteratively correcting the relative coordinates Ryckaert et al. (1977); Ciccotti et al. (1982), and, optionally, also the relative velocities Anderson (1983), using a relaxation procedure to ensure the constraints remain satisfied; the other tackles the problem by constructing a matrix representing the contributions of the constraints which, in effect, must be inverted at each timestep Edberg et al. (1986); Rapaport (1995), and which is subject to gradual drift requiring regular correction. Similar geometric constraints can be introduced to maintain constant bond angles as well, since it is often a reasonable approximation to assume that the angles between consecutive bonds along the backbone (or elsewhere) are unvarying. Such geometrical constraints have proved extremely useful, given the nature of the excitations present in the system: fluctuations in bond lengths, and sometimes also angles, tend to be of relatively small amplitude and high frequency, so that freezing them out of the dynamics permits a substantial increase in the allowed integration timestep. The amount of additional processing required for the constraints depends on their number $`n_c`$; the dependence is typically $`O(n_c)`$ for the iterative approach, but for the matrix approach it is $`O(n_c^3)`$, making the latter unsuitable for large problems.

If bonds lengths and angles are fixed, the only remaining internal DOFs are the dihedral angles, each defined in terms of a rotation about an axis lying along a bond, and affecting the relative orientation of the pair of bonds on either side. For reasons shrouded in history, dealing with this problem has been perceived as difficult, as indeed it is, if the problem is not addressed in a suitable manner. A significant advance in the methodology for dealing with dynamical problems involving internal coordinates occurred some years ago in the robotics field Rodriguez and Kreutz-Delgado (1992); Jain (1991), but with only the occasional exception, e.g., Bertsch et al. (1998), it appears to have gone unappreciated by the polymer simulation community at large. Because of the importance of this technique, the goal of which is to deal directly and economically with the internal DOFs, and since there is no reason why it should not be capable of replacing the various constraint-based approaches for most applications, a detailed treatment of the underlying theory is included in the paper.

This approach to the dynamics of linked bodies also requires solving the dynamics of individual rigid bodies. An alternative, recently described means Dullweber et al. (1997) of numerically dealing with the rigid-body equations of motion is discussed briefly; the method is based on rotation matrices, rather than on quaternions (or even Euler angles) that are generally used. The present formulation differs slightly from the original in regard to the reference frame in which the computations are carried out. The use of rotation matrices offers improved numerical stability, and since the method belongs to the leapfrog family of integrators, it means that simple leapfrog integration techniques can be used for the entire set of dynamical equations appearing in the problem.

## II Linked-body dynamics

### II.1 Chain coordinates

Consider a linear polymer chain whose monomers are joined by rigid bonds. In the discussion that follows, the terms “monomer”, “atom”, “site” and “joint” will be used interchangeably, as appropriate to the context, likewise “link” and “bond”. Bond lengths and angles are constant. If each torsional DOF is regarded as a mechanical joint associated with the site at one end of the link, with just a single rotational DOF, then the system is analogous to a basic problem in the field of robotic manipulators Rodriguez and Kreutz-Delgado (1992); Jain (1991).

The chain configuration is defined by the site positions $`\{𝒓_k\}`$, and if the bond vectors between adjacent sites are $`\{𝒃_k\}`$ then $`𝒓_{k+1}=𝒓_k+𝒃_k`$. The internal configuration of the chain can be specified by a set of bond rotation matrices $`\{𝑹_k\}`$. The transformation between the local coordinate frames attached to bonds $`k1`$ and $`k`$ ($`k1`$) involves a rotation through the bond angle $`\alpha _k`$ about the axis $`\widehat{𝒙}_{k1}`$, where $`\mathrm{cos}\alpha _k=\widehat{𝒃}_{k1}\widehat{𝒃}_k`$, followed by a rotation through the dihedral angle $`\theta _k`$ about the joint axis $`\widehat{𝒛}_{k1}`$. The matrix (actually its transpose) corresponding to this rotation is

$$𝑹_{k1,k}^T=\left(\begin{array}{ccc}\mathrm{cos}\theta _k& \mathrm{sin}\theta _k\mathrm{cos}\alpha _k& \mathrm{sin}\theta _k\mathrm{sin}\alpha _k\\ \mathrm{sin}\theta _k& \mathrm{cos}\theta _k\mathrm{cos}\alpha _k& \mathrm{cos}\theta _k\mathrm{sin}\alpha _k\\ 0& \mathrm{sin}\alpha _k& \mathrm{cos}\alpha _k\end{array}\right),$$
(1)
so that

$$𝑹_k^T=𝑹_0^T𝑹_{0,1}^T\mathrm{}𝑹_{k1,k}^T,$$
(2)
where $`𝑹_0^T`$ represents the orientation of the initial site and bond, and

$$𝒓_{k+1}=𝒓_k+|𝒃_k|𝑹_k^T\widehat{𝒛}.$$
(3)
In the present case, $`\{|𝒃_k|\}`$ and $`\{\alpha _k\}`$ are all constant, so that the only internal DOFs are those associated with $`\{\theta _k\}`$. Define $`\widehat{𝒉}_k`$ to be the rotation axis of the joint between bonds $`k1`$ and $`k`$ that is fixed in the frame of bond $`k1`$; in the present case $`\widehat{𝒉}_k\widehat{𝒛}_{k1}`$. Insofar as indexing is concerned, there are $`n_r`$ internal rotational joints (with labels $`1,\mathrm{},n_r`$), $`n_b=n_r+1`$ bonds ($`0,\mathrm{},n_r`$) and $`n_r+2`$ sites ($`0,\mathrm{},n_r+1`$). In order to completely specify the chain configuration, an additional joint is attached to the $`k=0`$ site, with three translational and three rotational DOFs (conceptually equivalent to a telescopic ball-and-socket joint); this joint is included in the formalism but will, eventually, be treated separately.

### II.2 Kinematic and dynamic relations

If $`𝒗_k`$ and $`𝝎_k`$ are the linear and angular velocities of site $`k`$, then the velocities and accelerations of adjacent sites are related by

$`𝝎_k`$ $`=`$ $`𝝎_{k1}+\widehat{𝒉}_k\dot{\theta }_k`$ (4)
$`𝒗_k`$ $`=`$ $`𝒗_{k1}+𝝎_{k1}\times 𝒃_{k1}`$ (5)
$`\dot{𝝎}_k`$ $`=`$ $`\dot{𝝎}_{k1}+\widehat{𝒉}_k\ddot{\theta }_k+𝝎_{k1}\times \widehat{𝒉}_k\dot{\theta }_k`$ (6)
$`\dot{𝒗}_k`$ $`=`$ $`\dot{𝒗}_{k1}+\dot{𝝎}_{k1}\times 𝒃_{k1}`$ (7)
$`+𝝎_{k1}\times (𝝎_{k1}\times 𝒃_{k1})`$
where $`1kn_r`$. While the mass elements of the chain are normally identified with the sites, here it is helpful to associate them with the bonds; if $`𝒓_k+𝒄_k`$ is the location of the center-of-mass of the atoms attached to bond $`k`$, then the center-of-mass acceleration of the bond is

$$\dot{𝒗_k^c}=\dot{𝒗}_k+\dot{𝝎}_k\times 𝒄_k+𝝎_k\times (𝝎_k\times 𝒄_k).$$
(8)
If $`𝒇_k`$ and $`𝒏_k`$ are the force and torque acting on bond $`k`$ across joint $`k`$, then the equations of motion are

$`𝓘_k\dot{𝝎}_k+𝝎_k\times (𝓘_k𝝎_k)`$ $`=`$ $`𝒏_k𝒏_{k+1}𝒄_k\times 𝒇_k`$ (9)
$`(𝒃_k𝒄_k)\times 𝒇_{k+1}+𝒏_k^e`$
$`m_k\dot{𝒗_k^c}`$ $`=`$ $`𝒇_k𝒇_{k+1}+𝒇_k^e,`$ (10)
where $`𝒇_k^e`$ and $`𝒏_k^e`$ are the externally applied force and torque; $`m_k`$ and $`𝓘_k`$ are the mass and moment of inertia of (the atoms associated with) the bond, the latter expressed in a space-fixed frame and relative to the center of mass of the bond. It is often convenient when dealing with rigid bodies to work in a center-of-mass frame Goldstein (1980); this is not the case here, and all vector components are expressed in the space-fixed coordinate frame. Rearrange the terms of Eqs. (9) and (10) to obtain relations between torques and forces on adjacent bonds,

$`𝒏_k`$ $`=`$ $`𝒏_{k+1}+𝒃_k\times 𝒇_{k+1}+m_k𝒄_k\times \dot{𝒗_k^c}`$ (11)
$`+𝓘_k\dot{𝝎}_k+𝝎_k\times (𝓘_k𝝎_k)𝒏_k^e𝒄_k\times 𝒇_k^e`$
$`𝒇_k`$ $`=`$ $`𝒇_{k+1}+m_k\dot{𝒗_k^c}𝒇_k^e,`$ (12)
and define the torque

$$t_k=\widehat{𝒉}_k𝒏_k$$
(13)
that acts along the axis $`\widehat{𝒉}_k`$ at joint $`k`$ and corresponds to the torsional interaction due to a twist around bond $`k1`$.

### II.3 Spatial operator formulation

Equations (4)–(7) can be expressed more concisely in terms of 6-component “spatial” vectors that combine the translational and rotational quantities. It is also convenient to represent certain vectors by means of antisymmetric matrices of form

$$\stackrel{~}{𝒖}=\left(\begin{array}{ccc}0& u_z& u_y\\ u_z& 0& u_x\\ u_y& u_x& 0\end{array}\right),$$
(14)
so that $`\stackrel{~}{𝒖}𝒗𝒖\times 𝒗`$. The resulting equations are

$`\left(\begin{array}{c}𝝎_k\\ 𝒗_k\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{I}& 0\\ \stackrel{~}{𝒃}_{k1}& \mathrm{I}\end{array}\right)\left(\begin{array}{c}𝝎_{k1}\\ 𝒗_{k1}\end{array}\right)+\left(\begin{array}{c}\widehat{𝒉}_k\\ 0\end{array}\right)\dot{\theta }_k`$ (23)
$`\left(\begin{array}{c}\dot{𝝎}_k\\ \dot{𝒗}_k\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{I}& 0\\ \stackrel{~}{𝒃}_{k1}& \mathrm{I}\end{array}\right)\left(\begin{array}{c}\dot{𝝎}_{k1}\\ \dot{𝒗}_{k1}\end{array}\right)+\left(\begin{array}{c}\widehat{𝒉}_k\\ 0\end{array}\right)\ddot{\theta }_k`$ (35)
$`+\left(\begin{array}{c}𝝎_{k1}\times \widehat{𝒉}_k\dot{\theta }_k\\ 𝝎_{k1}\times (𝝎_{k1}\times 𝒃_{k1})\end{array}\right)`$
or, equivalently,

$`\mathrm{V}_k`$ $`=`$ $`\varphi _{k1,k}^T\mathrm{V}_{k1}+\mathrm{H}_k^T\dot{\mathrm{W}}_k`$ (36)
$`\mathrm{A}_k`$ $`=`$ $`\varphi _{k1,k}^T\mathrm{A}_{k1}+\mathrm{H}_k^T\ddot{\mathrm{W}}_k+\mathrm{X}_k,`$ (37)
where $`\mathrm{V}_k`$ and $`\mathrm{A}_k`$ are examples of spatial vectors, and

$$\varphi _{k1,k}^T=\left(\begin{array}{cc}\mathrm{I}& 0\\ \stackrel{~}{𝒃}_{k1}& \mathrm{I}\end{array}\right).$$
(38)
The $`6\times 6`$ matrices $`\varphi _{k1,k}^T`$ and $`\varphi _{k,k+1}`$ (later) appear throughout the derivation, and their role is to propagate kinematic and dynamic information between joints. Several other new variables have been used:

$$\mathrm{H}_k^T=\left(\begin{array}{c}\widehat{𝒉}_k\\ 0\end{array}\right)$$
(39)
is a 6-component joint axis vector (in the more general case of a joint with $`d`$ DOFs, which the formalism is capable of handling, $`\mathrm{H}_k^T`$ would become a $`6\times d`$ matrix),

$$\mathrm{X}_k=\left(\begin{array}{cc}\stackrel{~}{𝝎}_{k1}& 0\\ 0& \stackrel{~}{𝝎}_{k1}\end{array}\right)\left(\begin{array}{c}\widehat{𝒉}_k\dot{\theta }_k\\ 𝒗_k𝒗_{k1}\end{array}\right)$$
(40)
is a 6-component spatial vector containing the remaining acceleration terms of the current site, and $`\dot{\mathrm{W}}_k\dot{\theta }_k`$. When used in vectors and matrices, $`\mathrm{I}`$ and $`0`$ denote unit and zero block submatrices of the implied size. The 6-component vectors, and most of the associated matrices, are shown in block capitals (to retain some similarity with Jain et al. (1993), $`\varphi `$, $`\psi `$ and $``$ are also used); no other special notation is needed since the variable types will be obvious from the context.

In a similar way, Eqs. (11)–(13) can be rewritten as

$`\left(\begin{array}{c}𝒏_k\\ 𝒇_k\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{I}& \stackrel{~}{𝒃}_k\\ 0& \mathrm{I}\end{array}\right)\left(\begin{array}{c}𝒏_{k+1}\\ 𝒇_{k+1}\end{array}\right)`$ (53)
$`+\left(\begin{array}{c}m_k𝒄_k\times \dot{𝒗_k^c}+𝓘_k\dot{𝝎}_k+𝝎_k\times (𝓘_k𝝎_k)\\ m_k\dot{𝒗_k^c}\end{array}\right)`$
$`\left(\begin{array}{c}𝒏_k^e+𝒄_k\times 𝒇_k^e\\ 𝒇_k^e\end{array}\right)`$
$`\left(\begin{array}{c}t_k\\ 0\end{array}\right)`$ $`=`$ $`\left(\begin{array}{c}\widehat{𝒉}_k\\ 0\end{array}\right)^T\left(\begin{array}{c}𝒏_k\\ 𝒇_k\end{array}\right)`$ (60)
or, equivalently,

$`\mathrm{F}_k`$ $`=`$ $`\varphi _{k,k+1}\mathrm{F}_{k+1}+\mathrm{M}_k\mathrm{A}_k+\mathrm{Y}_k`$ (61)
$`\mathrm{T}_k`$ $`=`$ $`\mathrm{H}_k\mathrm{F}_k.`$ (62)
Here Eq. (8) has been used, and

$$\mathrm{M}_k=\left(\begin{array}{cc}𝓘_km_k\stackrel{~}{𝒄}_k\stackrel{~}{𝒄}_k& m_k\stackrel{~}{𝒄}_k\\ m_k\stackrel{~}{𝒄}_k& m_k\mathrm{I}\end{array}\right),$$
(63)
is the symmetric, $`6\times 6`$ mass matrix; the 6-component vector

$$\mathrm{Y}_k=\left(\begin{array}{c}\stackrel{~}{𝝎}_k(𝓘_km_k\stackrel{~}{𝒄}_k\stackrel{~}{𝒄}_k)𝝎_k\\ m_k\stackrel{~}{𝝎}_k\stackrel{~}{𝝎}_k𝒄_k\end{array}\right)\left(\begin{array}{c}𝒏_k^e+𝒄_k\times 𝒇_k^e\\ 𝒇_k^e\end{array}\right)$$
(64)
contains the remaining force contributions. The identity $`𝒄_k\times (𝝎_k\times (𝝎_k\times 𝒄_k))=𝝎_k\times (𝒄_k\times (𝒄_k\times 𝝎_k))`$ was used in obtaining these expressions. In order to use the recurrence relations for $`\mathrm{V}_k`$, $`\mathrm{A}_k`$ and $`\mathrm{F}_k`$, the velocity and acceleration of the initial site, $`\mathrm{V}_0`$ and $`\mathrm{A}_0`$, must be provided, while the force associated with the site at the end of the final bond, $`\mathrm{F}_{n_r+1}`$, is zero, since there is no joint associated with that site.

The purpose of the recurrence relations in Eqs. (37) and (61) is to provide expressions for $`\{\ddot{\mathrm{W}}_k\}`$, which, together with $`\mathrm{A}_0`$, and assuming all the forces acting on the sites are known, can be integrated to solve for the chain dynamics; this is actually the opposite of the typical robotics problem, in which the goal is to determine the forces required to produce a particular robot arm trajectory.

### II.4 Stacked operators

Equations (36), (37), (61) and (62) can be rewritten in condensed, “stacked” form

$`\mathrm{V}`$ $`=`$ $`\varphi ^T\mathrm{V}+\mathrm{H}^T\dot{\mathrm{W}}`$ (65)
$`\mathrm{A}`$ $`=`$ $`\varphi ^T\mathrm{A}+\mathrm{H}^T\ddot{\mathrm{W}}+\mathrm{X}`$ (66)
$`\mathrm{F}`$ $`=`$ $`\varphi \mathrm{F}+\mathrm{MA}+\mathrm{Y}`$ (67)
$`\mathrm{T}`$ $`=`$ $`\mathrm{HF}`$ (68)
that combines the entire set of $`k`$ values. A quantity such as $`\mathrm{V}`$ containing all the $`\mathrm{V}_k`$ values for the chain is also referred to as a spatial vector, while, for example, the block matrix $`\varphi `$ containing all the $`\varphi _{k,k+1}`$ matrices is a spatial operator. The stacked formalism leads to a concise and elegant formulation of the problem, free from inundation by indices as is often the case in the robotics literature, e.g., Fu et al. (1987).

The spatial operator approach was originally developed for the case of a fixed initial bond Rodriguez and Kreutz-Delgado (1992) – the base in the example of a robot arm – for which $`\mathrm{V}_0=0`$, so that $`\dot{\mathrm{W}}=\mathrm{col}(\dot{\theta }_1,\mathrm{},\dot{\theta }_{n_r})`$ is a vector with just $`n_r`$ components, and the other vectors and matrices are sized accordingly. In order to remove the fixed-base restriction Jain and Rodriguez (1995), six extra DOFs are added to the problem by redefining $`\dot{\mathrm{W}}=\mathrm{col}(\mathrm{V}_0,\dot{\theta }_1,\mathrm{},\dot{\theta }_{n_r})`$ as a vector with $`n_r+6`$ components; likewise for $`\ddot{\mathrm{W}}`$. The size of the original $`6n_r\times n_r`$ block-diagonal matrix $`\mathrm{H}=\mathrm{diag}(\mathrm{H}_1,\mathrm{},\mathrm{H}_{n_r})`$ is increased to $`6(n_r+1)\times (n_r+6)`$ by including an extra $`6\times 6`$ block $`\mathrm{H}_0=\mathrm{I}`$, so that now $`\mathrm{H}=\mathrm{diag}(\mathrm{I},\mathrm{H}_1,\mathrm{},\mathrm{H}_{n_r})`$. The block-diagonal matrix $`\mathrm{M}`$ is of size $`6(n_r+1)\times 6(n_r+1)`$; $`\varphi `$ has the same size, and its only non-zero blocks are those to the immediate right of the diagonal, namely $`\{\varphi _{01},\mathrm{},\varphi _{n_r1,n_r}\}`$. Vectors $`\mathrm{V}`$, $`\mathrm{A}`$, $`\mathrm{F}`$, $`\mathrm{X}`$, and $`\mathrm{Y}`$, all have $`6(n_r+1)`$ components, e.g., $`\mathrm{V}=\mathrm{col}(\mathrm{V}_0,\mathrm{},\mathrm{V}_{n_r})`$, and $`\mathrm{T}`$ is organized in the same way as $`\dot{\mathrm{W}}`$, with $`n_r+6`$ components; $`\mathrm{T}_0=0`$ because the special $`k=0`$ joint exerts no torque. (Note that index order has been reversed from the original to make it more suitable for polymer use, and, for convenience, other aspects of the notation have been altered or simplified.)

The next step is to define the matrix

$$\mathrm{\Phi }=(\mathrm{I}\varphi )^1,$$
(69)
which is also used in the alternative form, $`\mathrm{\Phi }=\mathrm{\Phi }\varphi +\mathrm{I}`$; because $`\varphi ^{n_r+1}=0`$, Eq. (69) is equivalent to $`\mathrm{\Phi }=\mathrm{I}+\varphi +\varphi ^2+\mathrm{}+\varphi ^{n_r}`$, which is an upper-triangular block matrix whose elements, each a $`6\times 6`$ matrix, are

$$\mathrm{\Phi }_{ij}=\{\begin{array}{cc}\mathrm{I}\hfill & j=i\hfill \\ \varphi _{i,i+1}\hfill & j=i+1\hfill \\ \varphi _{i,i+1}\mathrm{}\varphi _{j1,j}\hfill & j>i+1.\hfill \end{array}$$
(70)
Then, in terms of $`\mathrm{\Phi }`$, Eqs. (65)–(68) reduce to

$`\mathrm{V}`$ $`=`$ $`\mathrm{\Phi }^T\mathrm{H}^T\dot{\mathrm{W}}`$ (71)
$`\mathrm{A}`$ $`=`$ $`\mathrm{\Phi }^T(\mathrm{H}^T\ddot{\mathrm{W}}+\mathrm{X})`$ (72)
$`\mathrm{T}`$ $`=`$ $`\ddot{\mathrm{W}}+\mathrm{H}\mathrm{\Phi }(\mathrm{M}\mathrm{\Phi }^T\mathrm{X}+\mathrm{Y}),`$ (73)
where

$$=\mathrm{H}\mathrm{\Phi }\mathrm{M}\mathrm{\Phi }^T\mathrm{H}^T.$$
(74)
While $`\mathrm{M}`$ is a sparse, $`6(n_r+1)\times 6(n_r+1)`$ block-diagonal matrix, $``$ is only of size $`(n_r+6)\times (n_r+6)`$, but, although it is typically much smaller, it is fully populated. In principle Eq. (73) can be numerically integrated to obtain $`\mathrm{W}`$, and this is one of the approaches actually used in solving the problem, but the computational effort required for evaluating $`^1`$ at each timestep to obtain $`\ddot{\mathrm{W}}`$ is of order $`O\left((n_r+6)^3\right)`$; for this reason such an approach is not practical for any but the shortest of chains. The alternative method, described below, requires a computational effort of order $`O(n_r)`$, together with what amounts to the inversion of a $`6\times 6`$ matrix; clearly this will prove to be a far more efficient approach, even for relatively small $`n_r`$.

### II.5 Inversion of the mass matrix

As a preliminary step in obtaining an explicit expression for $`^1`$ define Rodriguez and Kreutz-Delgado (1992) the $`6\times 6`$ matrix $`\mathrm{P}_k`$ in terms of $`\mathrm{M}_k`$ as

$$\mathrm{P}_k=\varphi _{k,k+1}(\mathrm{I}\mathrm{G}_{k+1}\mathrm{H}_{k+1})\mathrm{P}_{k+1}\varphi _{k,k+1}^T+\mathrm{M}_k.$$
(75)
In Eq. (75),

$`\mathrm{G}_k`$ $`=`$ $`\mathrm{P}_k\mathrm{H}_k^T\mathrm{D}_k^1`$ (76)
$`\mathrm{D}_k`$ $`=`$ $`\mathrm{H}_k\mathrm{P}_k\mathrm{H}_k^T,`$ (77)
where, for joints with a single DOF, $`\mathrm{G}_k`$ is a 6-component vector and $`\mathrm{D}_k`$ is a nonzero scalar; note also that $`\mathrm{P}_k`$ is symmetric. (The motivation for introducing $`\mathrm{P}_k`$ is explained in Rodriguez and Kreutz-Delgado (1992) and derives from the formal similarity of these equations with those used in the completely unrelated field of linear filtering.) Also define

$$\psi _{k,k+1}=\varphi _{k,k+1}(\mathrm{I}\mathrm{G}_{k+1}\mathrm{H}_{k+1})$$
(78)
and substitute this in Eq. (75). The stacked versions of Eqs. (75)–(78) are

$`\mathrm{P}`$ $`=`$ $`\psi \mathrm{P}\varphi ^T+\mathrm{M}`$ (79)
$`\mathrm{G}`$ $`=`$ $`\mathrm{PH}^T\mathrm{D}^1`$ (80)
$`\mathrm{D}`$ $`=`$ $`\mathrm{HPH}^T`$ (81)
$`\psi `$ $`=`$ $`\varphi (\mathrm{I}\mathrm{GH}).`$ (82)
Matrices $`\mathrm{P}`$ and $`\psi `$ are of size $`6(n_r+1)\times 6(n_r+1)`$, and $`\mathrm{G}`$ is $`(n_r+6)\times 6(n_r+1)`$ and block-diagonal (thus the product $`\mathrm{G}_{k+1}\mathrm{H}_{k+1}`$ is square). Matrix $`\mathrm{D}`$ is of size $`(n_r+6)\times (n_r+6)`$; its first $`6\times 6`$ diagonal block corresponds to $`\mathrm{D}_0`$, and the remaining $`n_r`$ diagonal elements are the scalars $`\mathrm{D}_k`$. From Eqs. (79) and (82),

$$\mathrm{M}=\mathrm{P}\varphi \mathrm{P}\varphi ^T+\varphi \mathrm{GHP}\varphi ^T,$$
(83)
and so, by using Eq. (69),

$$\mathrm{\Phi }\mathrm{M}\mathrm{\Phi }^T=\mathrm{P}+\mathrm{\Phi }\varphi \mathrm{P}+\mathrm{P}\varphi ^T\mathrm{\Phi }^T+\mathrm{\Phi }\varphi \mathrm{PH}^T\mathrm{D}^1\mathrm{HP}\varphi ^T\mathrm{\Phi }^T.$$
(84)
Substitute Eq. (84) in (74), then use $`\mathrm{GD}=\mathrm{PH}^T`$ from (80), together with (81), to obtain

$``$ $`=`$ $`\mathrm{HPH}^T+\mathrm{H}\mathrm{\Phi }\varphi \mathrm{PH}^T+\mathrm{HP}\varphi ^T\mathrm{\Phi }^T\mathrm{H}^T`$ (85)
$`+\mathrm{H}\mathrm{\Phi }\varphi \mathrm{PH}^T\mathrm{D}^1\mathrm{HP}\varphi ^T\mathrm{\Phi }^T\mathrm{H}^T`$
$`=`$ $`(\mathrm{I}+\mathrm{H}\mathrm{\Phi }\varphi \mathrm{G})\mathrm{D}(\mathrm{I}+\mathrm{H}\mathrm{\Phi }\varphi \mathrm{G})^T.`$
This alternative factorization of $``$ is a product of three $`(n_r+6)\times (n_r+6)`$ matrices, unlike Eq. (74) which involves nonsquare matrices.

It is now a straightforward matter to invert $``$. Use a special case of the Woodbury formula for the inverse of a matrix Press et al. (1992), $`(\mathrm{I}+\mathrm{Q}_1\mathrm{Q}_2)^1=\mathrm{I}\mathrm{Q}_1(\mathrm{I}+\mathrm{Q}_2\mathrm{Q}_1)^1\mathrm{Q}_2`$ to write

$$(\mathrm{I}+\mathrm{H}\mathrm{\Phi }\varphi \mathrm{G})^1=\mathrm{I}\mathrm{H}\mathrm{\Phi }(\mathrm{I}+\varphi \mathrm{GH}\mathrm{\Phi })^1\varphi \mathrm{G}.$$
(86)
By analogy with Eq. (69) for $`\mathrm{\Phi }`$, define $`\mathrm{\Psi }=(\mathrm{I}\psi )^1`$; then from Eqs. (82) and (69),

$$\mathrm{\Psi }^1=\mathrm{\Phi }^1+\varphi \mathrm{GH},$$
(87)
so that $`(\mathrm{I}+\mathrm{H}\mathrm{\Phi }\varphi \mathrm{G})^1=\mathrm{I}\mathrm{H}\mathrm{\Psi }\varphi \mathrm{G}`$. Thus the inverse of Eq. (85) is

$$^1=(\mathrm{I}\mathrm{H}\mathrm{\Psi }\varphi \mathrm{G})^T\mathrm{D}^1(\mathrm{I}\mathrm{H}\mathrm{\Psi }\varphi \mathrm{G}),$$
(88)
and so, from Eq. (73),

$`\ddot{\mathrm{W}}`$ $`=`$ $`(\mathrm{I}\mathrm{H}\mathrm{\Psi }\varphi \mathrm{G})^T\mathrm{D}^1(\mathrm{I}\mathrm{H}\mathrm{\Psi }\varphi \mathrm{G})`$ (89)
$`\times [\mathrm{T}\mathrm{H}\mathrm{\Phi }(\mathrm{M}\mathrm{\Phi }^T\mathrm{X}+\mathrm{Y})]`$
$`=`$ $`(\mathrm{I}\mathrm{H}\mathrm{\Psi }\varphi \mathrm{G})^T\mathrm{D}^1`$
$`\times [\mathrm{T}\mathrm{H}\mathrm{\Psi }(\varphi \mathrm{GT}+\mathrm{M}\mathrm{\Phi }^T\mathrm{X}+\mathrm{Y})],`$
where Eq. (87) is used in simplifying $`\mathrm{H}(\mathrm{I}\mathrm{\Psi }\varphi \mathrm{GH})\mathrm{\Phi }=\mathrm{H}\mathrm{\Psi }`$. To eliminate $`\mathrm{\Psi }`$, first rewrite Eq. (89) as

$$(\mathrm{I}+\mathrm{H}\mathrm{\Phi }\varphi \mathrm{G})^T\ddot{\mathrm{W}}=\mathrm{D}^1[\mathrm{T}\mathrm{H}\mathrm{\Psi }(\varphi \mathrm{GT}+\mathrm{M}\mathrm{\Phi }^T\mathrm{X}+\mathrm{Y})].$$
(90)
Next, use Eq. (79) with (69) to get

$`\mathrm{\Psi }\mathrm{M}\mathrm{\Phi }^T`$ $`=`$ $`\mathrm{\Psi }\mathrm{P}(\varphi ^T\mathrm{\Phi }^T+\mathrm{I})\mathrm{\Psi }\psi \mathrm{P}\varphi ^T\mathrm{\Phi }^T`$ (91)
$`=`$ $`\mathrm{\Psi }\mathrm{P}+\mathrm{P}\varphi ^T\mathrm{\Phi }^T.`$
Then, using the transpose of Eq. (80), it follows that

$$(\mathrm{I}+\mathrm{H}\mathrm{\Phi }\varphi \mathrm{G})^T\ddot{\mathrm{W}}=\mathrm{D}^1\mathrm{E}\mathrm{G}^T\varphi ^T\mathrm{\Phi }^T\mathrm{X},$$
(92)
in which the force-like quantities

$`\mathrm{E}`$ $`=`$ $`\mathrm{T}\mathrm{HZ}`$ (93)
$`\mathrm{Z}`$ $`=`$ $`\mathrm{\Psi }(\varphi \mathrm{GT}+\mathrm{PX}+\mathrm{Y})`$ (94)
have been defined. Rearranging Eq. (92) and using the expression for $`\mathrm{A}`$ given in Eq. (72) leads to

$`\ddot{\mathrm{W}}`$ $`=`$ $`\mathrm{D}^1\mathrm{E}\mathrm{G}^T\varphi ^T\mathrm{\Phi }^T(\mathrm{H}^T\ddot{\mathrm{W}}+\mathrm{X})`$ (95)
$`=`$ $`\mathrm{D}^1\mathrm{E}\mathrm{G}^T\varphi ^T\mathrm{A}.`$
It is also possible to eliminate $`\mathrm{\Psi }`$ from Eq. (94) by substituting $`\mathrm{T}`$ from (93) to get $`(\mathrm{I}\mathrm{\Psi }\varphi \mathrm{GH})\mathrm{Z}=\mathrm{\Psi }(\varphi \mathrm{GE}+\mathrm{PX}+\mathrm{Y})`$, and then using Eq. (87) to obtain

$$\mathrm{Z}=\mathrm{\Phi }(\varphi \mathrm{GE}+\mathrm{PX}+\mathrm{Y}).$$
(96)
Explicit forms for the new recurrence relations embodied in Eqs. (95) and (96) are obtained by using Eq. (69) and reintroducing the $`k`$ indices –

$`\mathrm{Z}_k`$ $`=`$ $`\varphi _{k,k+1}(\mathrm{Z}_{k+1}+\mathrm{G}_{k+1}\mathrm{E}_{k+1})+\mathrm{P}_k\mathrm{X}_k+\mathrm{Y}_k`$ (97)
$`\ddot{\mathrm{W}}_k`$ $`=`$ $`\mathrm{D}_k^1\mathrm{E}_k\mathrm{G}_k^T\varphi _{k1,k}^T\mathrm{A}_{k1}.`$ (98)
These recurrence relations are used in opposite $`k`$-directions; they succeed in providing the required results without the need for explicit evaluation of the matrix inverse $`^1`$ as implied by Eq. (73). It is for this reason that the method has not been referred to as an “inverse matrix method”, a term sometimes seen in the literature, but rather a “rigid link” method, a far more apt descriptor.

The expressions given here describe the entire chain, but, provided the end joints are handled correctly, these results can be used for linear segments that form part of a larger assembly, allowing more complicated tree-like structures to be treated. Furthermore, while the above formulation deals with the simplest case of a linear chain with a single torsional DOF per joint, it is readily extended to more complex joints, enabling, for example, the constant bond-angle condition to be eliminated by allowing two DOFs at each joint (an alternative would be to decompose an individual joint into two coincident joints each with a single DOF).

## III Simulation techniques

### III.1 Linked-chain equations of motion

The recurrence relations used to propagate velocities, forces, and accelerations along the chain are as follows: The (translational and rotational) velocities $`\mathrm{V}_k`$ are obtained by starting with $`\mathrm{V}_0`$ and iterating Eq. (36),

$$\mathrm{V}_k=\varphi _{k1,k}^T\mathrm{V}_{k1}+\mathrm{H}_k^T\dot{\mathrm{W}}_k,k=1,\mathrm{},n_r.$$
(99)
The forces (and torques), as represented by $`\mathrm{E}_k`$, together with the matrices $`\mathrm{D}_k`$ and $`\mathrm{G}_k`$, are obtained by iterating Eqs. (75) and (97). For computational convenience, new quantities $`\mathrm{A}_k^{}`$ and $`\mathrm{Z}_k^{}`$ are introduced; then, starting with $`\mathrm{P}_{n_r+1}=0`$ and $`\mathrm{Z}_{n_r+1}^{}=0`$,

$$\begin{array}{ccc}\mathrm{P}_k& =& \varphi _{k,k+1}(\mathrm{I}\mathrm{G}_{k+1}\mathrm{H}_{k+1})\hfill \\ & & \times \mathrm{P}_{k+1}\varphi _{k,k+1}^T+\mathrm{M}_k\hfill \\ \mathrm{D}_k& =& \mathrm{H}_k\mathrm{P}_k\mathrm{H}_k^T\hfill \\ \mathrm{G}_k& =& \mathrm{P}_k\mathrm{H}_k^T\mathrm{D}_k^1\hfill \\ \mathrm{Z}_k& =& \varphi _{k,k+1}\mathrm{Z}_{k+1}^{}+\mathrm{P}_k\mathrm{X}_k+\mathrm{Y}_k\hfill \\ \mathrm{E}_k& =& \mathrm{T}_k\mathrm{H}_k\mathrm{Z}_k\hfill \\ \mathrm{Z}_k^{}& =& \mathrm{Z}_k+\mathrm{G}_k\mathrm{E}_k\hfill \end{array}\}k=n_r,\mathrm{},0.$$

Finally, the values of $`\ddot{\mathrm{W}}_k`$ (or $`\ddot{\theta }_k`$) are determined by starting with $`\mathrm{A}_0`$ (its evaluation is discussed below), and iterating Eqs. (37) and (98),

$$\begin{array}{ccc}\mathrm{A}_k^{}& =& \varphi _{k1,k}^T\mathrm{A}_{k1}\hfill \\ \ddot{\mathrm{W}}_k& =& \mathrm{D}_k^1\mathrm{E}_k\mathrm{G}_k^T\mathrm{A}_k^{}\hfill \\ \mathrm{A}_k& =& \mathrm{A}_k^{}+\mathrm{H}_k^T\ddot{\mathrm{W}}_k+\mathrm{X}_k\hfill \end{array}\}k=1,\mathrm{},n_r.$$
(100)
These recurrence relations, which are readily transformed into a suitable computer program, imply a series of operations (multiplications and additions) involving $`6\times 6`$ matrices and 6-component vectors, but the total computational effort is only of order $`O(n_r)`$.

Recall that the $`k=0`$ joint has six DOFs, and also that $`\mathrm{H}_0=\mathrm{I}`$, $`\mathrm{X}_0=0`$, and $`\ddot{\mathrm{W}}_0=\mathrm{A}_0`$. Now, because $`\mathrm{A}_1=0`$, it follows from Eq. (100) that $`\mathrm{A}_0=\mathrm{D}_0^1\mathrm{E}_0`$, and since $`\mathrm{T}_0=0`$,

$$\mathrm{D}_0\mathrm{A}_0=\mathrm{Z}_0,$$
(101)
where both $`\mathrm{D}_0`$ and $`\mathrm{Z}_0`$ have already been determined (above). Thus $`\mathrm{A}_0`$ can be evaluated numerically by solving the set of six linear equations contained in Eq. (101) using the standard LU method Press et al. (1992); the computational effort required for this initial joint is fixed and independent of $`n_r`$.

### III.2 Leapfrog integration and rigid-body equations

The familiar leapfrog method for integrating the MD translational equations of motion – which is algebraically equivalent to the Verlet method Verlet (1967) – is usually expressed in a form where the coordinates and velocities are evaluated at alternate half-timesteps Rapaport (1995). This minor inconvenience can be avoided by using a slightly modified form that breaks the integration procedure for a single timestep into two parts: Prior to computing the latest acceleration ($`𝒂`$) values, update the velocities ($`𝒗`$) by a half timestep using the previous accelerations, and then update the coordinates ($`𝒓`$) by a full timestep using these intermediate velocity values,

$`𝒗(t+h/2)`$ $`=`$ $`𝒗(t)+(h/2)𝒂(t)`$ (102)
$`𝒓(t+h)`$ $`=`$ $`𝒓(t)+h𝒗(t+h/2).`$ (103)
In the case of the polymer chain, this procedure is applied to the translation coordinates of the $`k=0`$ site and (in scalar form) to each of the dihedral angles $`\theta _k`$; the treatment of the angular coordinates associated with the $`k=0`$ site, below, employs a related approach for dealing with the rotational equations. Next, use the new coordinates (and velocities if needed) to compute the latest acceleration values, then update the velocities over the second half timestep,

$$𝒗(t+h)=𝒗(t+h/2)+(h/2)𝒂(t+h).$$
(104)
In the linked-chain formulation, the initial bond of the chain is treated as a rigid body; the influence of the rest of the chain on it has already been taken into account and is contained in the force and torque transmitted through the first internal joint. There are a number of ways of describing the orientation of a rigid body Goldstein (1980): Euler angles have proved very useful for analytic purposes because of their intuitive nature, but owing to a potentially singular matrix that appears in the equations of motion they are not the preferred method for dealing with numerical problems. Quaternions have achieved popularity because of their singularity-free nature, but their normalization must be preserved against a small but persistent numerical drift Evans (1977); Rapaport (1995). A more recently proposed alternative is to regard the complete rotation matrix as the dynamical variable; this is the representation that will be used here, since the integration scheme Dullweber et al. (1997) – which is based on operator splitting and maintains time-reversibility – is just another instance of the leapfrog method.

In the original description Dullweber et al. (1997), vector components were expressed in the principal-axis frame of the body. Since the chain dynamical problem as a whole is solved in the space-fixed frame, the corresponding form of the rotational equations will be described here. If $`𝑹`$ denotes the rotation matrix of a rigid body, then the first part of the leapfrog integration step consists of a half-timestep update of the angular velocities,

$$𝝎(t+h/2)=𝝎(t)+(h/2)𝜶(t),$$
(105)
where $`𝜶\dot{𝝎}`$, followed by a full-timestep update of $`𝑹`$ using a symmetric product of matrices describing a series of small partial rotations,

$$𝑹^T(t+h)=𝑼_1𝑼_2𝑼_3𝑼_2𝑼_1𝑹^T(t),$$
(106)
where, for convenience, the transpose of $`𝑹`$ is treated. Note that for the linked chain, the rigid body is associated with the $`k=0`$ joint, so that $`𝑹𝑹_0`$. Each of the matrices

$$𝑼_1=𝑼_x(\omega _xh/2),𝑼_2=𝑼_y(\omega _yh/2),𝑼_3=𝑼_z(\omega _zh)$$
(107)
describes a rotation about a single axis and is evaluated in the space-fixed frame. For small angles, they can be approximated in a way that preserves orthogonality, e.g.,

$$𝑼_x(\theta )=\left(\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\theta & \mathrm{sin}\theta \\ 0& \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& \frac{1\theta ^2/4}{1+\theta ^2/4}& \frac{\theta }{1+\theta ^2/4}\\ 0& \frac{\theta }{1+\theta ^2/4}& \frac{1\theta ^2/4}{1+\theta ^2/4}\end{array}\right).$$
(108)
The second part of the leapfrog step is

$$𝝎(t+h)=𝝎(t+h/2)+(h/2)𝜶(t+h).$$
(109)
In the case of a single rigid body, the angular acceleration is determined from the torque $`𝝉`$, namely $`𝜶(t+h)=𝓘^1𝝉(t+h)`$, whereas for the linked chain this treatment is only required for the $`k=0`$ joint, and $`𝜶`$ is obtained by solving Eq. (101); the reason rigid bodies are usually treated in the body-fixed principal-axes frame is to ensure the diagonality of $`𝓘`$, a consideration that is not relevant here.

The complete procedure for a single timestep can be summarized as the following sequence of operations: integrate (first part) to obtain base velocities and coordinates, and joint angular velocities and angles; determine site velocities, Eq. (99); evaluate site coordinates, Eqs. (1)–(3); compute external forces and torques, and other necessary quantities; determine joint forces, Eq. (III.1); solve Eq. (101) for the base acceleration; determine joint accelerations, Eq. (100); integrate (second part) to obtain base velocities and joint angular velocities.

### III.3 Polymer chain model

Two kinds of interactions are required in this model – excluded volume and torsion. The former is provided by a pair interaction that prevents overlap of the atoms (or atom groups) located at the chain sites. Here a simple soft-sphere repulsion, based on the Lennard-Jones potential with a short-range cutoff, is all that is required: for a pair of atoms located at $`𝒓_i`$ and $`𝒓_j`$, where $`𝒓_{ij}=𝒓_i𝒓_j`$, and $`r_{ij}=|𝒓_{ij}|`$, the potential is

$$u_{ss}(r_{ij})=\{\begin{array}{cc}4ϵ[(r_{ij}/\sigma )^{12}(r_{ij}/\sigma )^6]& r_{ij}<r_c\\ 0& r_{ij}r_c\end{array}$$
(110)
with a cutoff $`r_c=2^{1/6}\sigma `$ (nearby pairs of atoms that are prevented from approaching too closely because of geometrical restrictions need not be considered). Should a pairwise attraction between particular pairs of distant chain atoms be required (as will be the case later on), it can be obtained from Eq. (110) by simply increasing $`r_c`$. The pair forces derived from this potential, and their associated torques, contribute to $`𝒇_k^e`$ and $`𝒏_k^e`$ in Eqs. (9) and (10).

The torsional potential associated with the dihedral angles $`\theta _k`$ has the simple form

$$u_t(\theta _k)=u_k\mathrm{cos}(\theta _k\theta _k^{(0)}),$$
(111)
where $`\theta _k^{(0)}`$ is the dihedral angle that produces a ground state having the correct helical twist, and $`u_k`$ is the interaction strength. The torque appearing in Eq. (13) is

$$t_k=u_k\mathrm{sin}(\theta _k\theta _k^{(0)}),$$
(112)
a result whose simplicity stands in sharp contrast to the intricate vector algebra associated with torque calculations when working in Cartesian coordinates Rapaport (1995).

For the chains considered here it is assumed all $`|𝒃_k|=b`$, $`\alpha _k=\alpha `$, $`\theta _k^{(0)}=\theta ^{(0)}`$ and, except for the later twin-helix studies where selected $`u_k=0`$, all $`u_k=u^{(0)}`$. Since the torsion also acts at the first internal joint, it is necessary to add an extra site and bond to the chain (effectively with an index “$`1`$”) to make this torsion term meaningful; the first three sites of the modified chain form a rigid unit (the extra bond does not alter the preceding analysis) and the chain length is increased by unity.

A spherical mass element (with a finite moment of inertia about its own center of mass) is associated with each site; for bonds with $`k>0`$, the mass is attached to the far ($`k+1`$ site) of the bond, while the $`k=0`$ bond, as explained above, has three masses associated with it. The components of the inertia tensor in Eqs. (63) and (64) are

$$(𝓘_k)_{ij}=\{\begin{array}{cc}_{\kappa k}m_\kappa (r_\kappa ^2r_{\kappa i}^2)\hfill & i=j\hfill \\ _{\kappa k}m_\kappa r_{\kappa i}r_{\kappa j}\hfill & ij\hfill \end{array},$$
(113)
where the sum (or volume integral) is over all mass elements $`\kappa `$ fixed to bond $`k`$, and coordinates are relative to the center of mass of each bond in the space-fixed frame.

### III.4 Order parameter

While the appearance of an ordered helical structure, even one with the occasional defect, is easily recognized visually, in order to facilitate statistical analysis of the behavior it is important to be able to quantify the degree of order present in the chain. Let $`𝒅_k=𝒃_{k1}\times 𝒃_k`$, then

$$S=\frac{1}{n_r}\underset{k=1}{\overset{n_r}{}}\widehat{𝒅}_k$$
(114)
defines an order parameter that measures the long-range order present in the folded structure based on the orientation of the helical turns; for a single, well-formed helix, $`S`$ should have a value close to unity. A slightly modified version of $`S`$ will be introduced later for studying twin-helix structures.

This definition of $`S`$ is particularly useful for detecting structures consisting of two or more helical domains with axes aligned in different directions due to a localized defect of the type seen in helically wound telephone and electrical cords. Since the correct helicity (or “handedness”) is built into the interactions, it is unlikely that segments of opposite helicity will independently nucleate at separate locations but, as the chain collapses, individual turns with the wrong twist can become trapped in the structure. These defects are capable of traveling along the chain, but this is a slow process, and the direction of motion is random unless close to the chain end. There are instances where the definition of $`S`$ in Eq. (114) can give an incomplete picture; if a wrong turn occurs very close to the chain end, its effect on $`S`$ will be minimal, and even a perfectly formed helix is subject to low frequency bending motion. Other order parameters can be defined that are of a more short-range nature; for example, a simple count of the number of pairs of chain sites lying within a specified range (i.e., the number of “contacts” between adjacent turns of the helix) divided by the maximum possible value, but for long chains the tolerance in the threshold required to accommodate thermal fluctuations might allow significant changes in the helical-axis direction to go undetected.

## IV Results

### IV.1 Simulation details

The most prominently recognizable structural motif found in proteins is the ($`\alpha `$-) helix. A helix, because of its uniformity along the longitudinal axis, is a particularly simple structure to specify, and both Monte Carlo and MD helix-folding simulations based on the complex potentials designed for protein modeling have been carried out, e.g., Rapaport and Scheraga (1981); Bertsch et al. (1998). Since the complexity of these potentials contributes little to a basic understanding of generic folding phenomena, the present simulations are based on the much simpler model and potentials described previously. Indeed, an analogous approach has been employed experimentally Nelson et al. (1997) in a study of helix formation in synthesized nonbiological chain molecules, where the interactions are simpler than in proteins (in particular, there are no hydrogen bonds).

The importance of examining simple structures, such as the helix, is that the process by which ordered arrangements emerge from randomly coiled states is likely to capture something of the essence of real protein folding, such as the cooperativity of the folding process, the role of nucleation sites, the degree to which folding is able to proceed to completion, and the steric and topological effects of excluded volume. The key question, of course, is whether the intrinsic timescales of these processes are sufficiently small for simulation to be computationally feasible, an issue addressed by the results presented here. While the present model is admittedly a mere caricature of the detailed models normally employed in studies of individual proteins, it has two undeniable advantages, namely a known native ground state compatible with the interactions, and sufficiently modest computational requirements that MD simulation is able to encompass the time interval required for major conformational change. More complex protein structures also display certain common characteristics, and ought to be accessible to simulations of this type; it is, however, essential to eliminate any ambiguity from the ground state, something that nature itself has presumably achieved in the interests of efficiency and reliability.

Each simulation run considers a single chain constructed as described earlier. The absence of a solvent, apart from changing the timescales, should not alter the outcome; indeed many, if not most, protein simulations avoid introducing an explicit solvent for reasons of computational efficiency. The simulation is begun at a relatively high temperature, so that the kinetic energy is sufficiently large to surmount the torsional potential barriers. The initial chain configuration is a large loop extending across the simulation cell, with a very slight helicity to prevent any overlap; initial dihedral angles are chosen so that locally, the conformation is almost a planar zigzag state. The joint angular velocities are assigned random values corresponding to the starting temperature, and memory of this initial state rapidly vanishes early in the simulation. The temperature is gradually reduced by a factor slightly less than unity at regular intervals until, towards the end of the run, very little kinetic energy remains in the system. The simulation region is bounded by hard, reflecting walls; while there are occasional wall collisions, this has little influence on the overall behavior. (The alternative would be to use periodic boundaries, which for a simulation cell not large enough to contain the chain in a fully-stretched state, would be subject to chain wraparound effects; while these are also unlikely to affect the overall behavior, they can prove visually confusing given the importance of computer-generated visualization in this work.)

The gradual cooling that is imposed throughout the run plays several distinct roles. During the early stage it is used to drive the chain from a totally random state to one in which the torsional potential begins to have some influence over the dihedral angles. Then, as the temperature is reduced further, an increasing degree of local order emerges and precursors to long-range order appear, either as a consequence of the merging of separate ordered domains, or the spread of order from a nucleation region (or a combination of both processes); during this stage the imposed cooling performs a task normally the responsibility of the solvent, namely the removal of excess potential energy as the chain evolves towards states of lower energy. Once the chain has reached a state consisting mainly of helical segments, possibly separated by small misfolded regions that have become trapped, the purpose of further temperature reduction is to gradually freeze out thermal fluctuations – without further major structural change – in order to allow evaluation of the long-range order parameter $`S`$ (the measure of success of the folding process); the latter part of this cooling stage is not intended to imitate any real physical process.

The simulations use standard, reduced MD units, in which all distances and energies are expressed in terms of the Lennard-Jones parameters $`\sigma `$ and $`ϵ`$ respectively; mass is expressed in terms of the monomer mass $`m`$ and, consequently, the unit of time is $`\sqrt{m\sigma ^2/ϵ}`$. Temperature and energy are made numerically identical by setting the Boltzmann constant $`k_\mathrm{B}`$ to unity. In terms of these units the parameters used in the runs are as follows: The bond length $`b=1.3`$, a value sufficiently short to prevent the chain crossing itself, the bond angle $`\alpha `$ and the preferred dihedral angle $`\theta ^{(0)}`$ are chosen to produce helices with periodicity six, and the torsional potential strength $`u^{(0)}=5`$. In the studies of twin helices, the cutoff in the attractive interaction, based on Eq. (110), occurs at $`r_c=2.2`$. The initial temperature is 4 (corresponding to a kinetic energy per DOF of 2) and the final temperature is $`10^3`$; temperature is reduced by rescaling all velocities and angular velocities by a factor $`f_T`$ every 4000 timesteps, with $`f_T=0.95`$ or 0.97. The runs reported here are each of length 4–8 $`\times 10^5`$ steps; the integration timestep (in MD units) is $`h=4\times 10^3`$. In order to produce reliable statistics, a large number of runs were carried out for each case studied; the runs differed in the choice of initial random values for $`\{\dot{\theta }_k\}`$.

### IV.2 Folding to a single helix

Measurements were made of the long-range order parameter $`S`$ and the total energy, the latter a sum over contributions from the soft-sphere pair interactions, Eq. (110), the torsional terms, Eq. (111), and the kinetic energy. The measurements involved 400 independent runs for each of several chain lengths $`L`$ and different cooling rates. These quantitative results were complemented by an interactive graphical version of the simulation program that provided real-time visual monitoring of the folding process; in addition to learning about any potential obstructions to complete folding, the ability to observe chains directly also helped when choosing a cooling rate sufficiently fast for folding to proceed to completion, but not too fast for an excessive number of defects to become trapped in the nascent structures.

The viability of the underlying approach depends on whether it can actually produce correctly structured helices. The first series of results measures the fraction of chains that successfully fold into a helical state, and the manner in which the success rate depends on $`L`$ and the cooling rate. A summary appears in Table 1; $`L`$ ranges from 18 to 90, which, since the helix period is six, corresponds to 3–15 full helical turns.

Owing to the large number of runs it is not possible to provide a detailed history of each, so a quantitative measure of folding success must be introduced. A successfully folded helix is deemed to be one for which $`S>0.88`$ at least once during the last $`1.2\times 10^5`$ steps of the run (measurements are made every 4000 steps); by this stage of the run the system has reached a comparatively low temperature, so that further substantial conformational changes are unlikely. Visual analysis confirms that, for the cases considered, this threshold for $`S`$ provides a quite reliable estimator; it tends to be sensitive to defects in the helical structure, while allowing for the fact that a properly folded helix may still has some residual curvature along its major axis.

It is clear from Table 1 that a high success rate for helix production is achieved. Two trends are apparent in the results, neither of them unexpected: for a given $`f_T`$, longer chains are less likely to fold properly than shorter chains, and, for a given $`L`$, a larger $`f_T`$ (corresponding to slower cooling) raises the success rate. Thus the longer the chain, the slower the desired cooling rate; additional runs with faster cooling confirm this observation. The longest of the chains folds to a helix with 15 turns, which, considering the potential for defects, represents a significant victory of energy over entropy.

The rate at which chains approach the helically ordered state can be studied by monitoring the mean values of $`S`$, as well as the negative of the total energy (which is dominated by the torsional component when in the folded state); these quantities provide measures of the long- and short-range order, respectively. The results, normalized per DOF, for 54- and 90-site chains, averaged over all 400 runs, are shown in Figs. 3 and 4. The overall results are divided into two groups, depending on whether the chain is classified as having folded successfully or not, and error bars indicate the standard deviations of the measurements. In each case it is the upper curve that represents the average for the successfully folded chains, and it has the smaller error bars.

An alternative estimate of the rate at which folding proceeds is based on the time dependence of the fraction of chains in a helical state, relative to all those that eventually succeed in reaching this state. This provides information about when, assuming a chain folds successfully, the appearance of helical order actually occurs. Fig. 5 shows these cumulative distributions for the different $`L`$, each at the slowest cooling rate considered. The cooling rate clearly affects the results, as can be seen from the separation of the two groups of curves that are based on different rates (see Table 1); for a given cooling rate, the folding speed tends to drop as the chains become longer (the slight crossover of the $`L=72`$ and 90 curves is probably not significant).

A more detailed examination of final-state conformations is based on histograms of the $`S`$-distribution, using measurements made over the last $`1.2\times 10^5`$ steps. These results appear in Fig. 6 for several $`L`$ values (at the slowest cooling rate). There are two separate curves for each $`L`$, one showing the spread of $`S`$ for those chains that satisfied the folding criterion at least once during this measurement period, and a broader, much lower curve for the chains that did not. The former set of distributions become broader with increasing $`L`$; there are several contributing causes for this, including slower folding rates (all the runs were of equal length), chains not managing to fold successfully but having one or more intermediate $`S`$ values which passed the test, and the increasing effect of bending along the helical axis. The latter, broader distributions (scaled up by a factor of 10 to make the details visible) are due both to the many defective structures possible, each with its own spread of $`S`$ values, and also to the defects themselves reducing the structural rigidity and so raising the susceptibility to slow thermal vibration. Only for the longest chains is there any overlap of the curves, and even then it is minimal.

The best way to follow the folding process is by viewing animated sequences of images taken at various points during the run; some sequences can actually be generated while running the simulation interactively, if the computations proceed sufficiently rapidly. Here, due to the limitations of the printed page, a selection of static images must suffice. One could attempt a verbal description of what transpires but, as was the case in Bertsch et al. (1998), there are no obvious features shared by the individual folding trajectories. Even if certain common characteristics do exist, the strong random conformational fluctuations make their observation difficult; a systematic, quantitative means for identifying pathways, that extends ideas used for equilibrium states Crippen and Ohkubo (1998), might prove helpful in this task.

Figure 1, which appeared early in the paper, shows an image of a typical, well-formed, almost straight helix obtained in one of the runs, while Fig. 2 shows a random chain configuration observed near the start of a run, both for $`L=90`$ chains. For clarity, these and subsequent pictures represent the chains by their tubular envelopes, rather than by showing individual, partially overlapped spheres representing the monomers. The tube thickness corresponds to unit diameter, slightly less than the soft-sphere interaction cutoff to ensure that a small amount of space remains visible between adjacent turns of the helix. The maximum amount of residual curvature that was observed in the backbones of properly folded helices is apparent from Fig. 7; this example actually meets the criterion for folding success (as indeed it should).

While the majority of runs (85% for chains with $`L=90`$, and an even higher proportion for smaller $`L`$, see Table 1) result in a correctly folded helix, examination of the kinds of defects that appear in the final states of those runs that fail to fold properly is an informative exercise. The first such picture, Fig. 8, is of an $`L=90`$ chain with two helical regions separated by a single defect. The defect is essentially a single loop of the helix with a reverse fold that became frozen in place during the cooling process. This is the most frequent type of defect, and its location can be anywhere in the chain, even right at the end.

Less frequent are chains with two spatially separated defects, as shown in Fig. 9. More extreme, but very rare examples of other kinds of defects appear in Figs. 10 and 11. These show what can happen when the chain starts to become entangled with itself; in the first case the problem is localized, but in the second (the only example of its kind observed) there is a relatively large loop trapped by the entanglement. The fact that these defects are relatively infrequent, and that even this low level of failure can be reduced by lowering the cooling rate still further, attest to the robustness and reliability of the folding process.

### IV.3 Folding to a pair of helices

The helix formation described above is obtained in a study of homopolymers where, due to the uniformity of the chain, the only kind of repeating structure that can be produced is the helix (the planar zigzag conformation is a degenerate case). In a protein context this kind of structure is classed as secondary because helices often serve as structural components in more complex assemblies. Globular proteins are characterized by at least one additional level in the structural hierarchy, namely tertiary structure. A model capable of demonstrating tertiary structure requires at least some differentiation among the chain sites that breaks the translational invariance. Building on the helix-forming model investigated here, the next stage of complexity is a packed assembly of helices, a structure that incorporates both the secondary and tertiary levels of the hierarchy. The simplest way to design such a structure is to include non-local, attractive forces between selected pairs of chain sites; here, the pairs involved are located at chain positions that will be brought into proximity following the collapse into a state with two adjacent helices aligned in antiparallel directions. Such highly specific interactions are reminiscent of an approach used for lattice protein models Ueda et al. (1978); the overall simplicity should be contrasted with the highly detailed model, complete with solvent, used in an MD study of the unfolding of a three-helix bundle Boczko and Brooks III (1995).

Choosing the interactions to produce a twin-helix structure is accomplished as follows: For a homogeneous chain with $`L=n_b+1`$ sites, in which the periodicity of the helix is $`p`$, the ground state consists of $`n_t=L/p`$ turns. Now assume that $`n_t`$ is an odd number and choose the interactions appropriate for a pair of adjacent helices, each with $`(n_t1)/2`$ turns, joined by a “bridging” chain segment of length $`p`$. All that remains is to identify the pairs of sites in the two helical regions that must attract; these are just neighboring sites in adjacent turns of one of the helical segments, matched with the corresponding sites, in reverse order, of the other. The strength of the attractive potential responsible for the tertiary structure, which is based on Eq. (110), is $`0.2u^{(0)}`$; it is weaker than the torsion, but the question of whether tertiary structure formation is the beneficiary or the cause of secondary structure formation Dill et al. (1995) is not addressed here. In these exploratory computations, the torsional interactions along the bonds in the bridging segment are set to zero for simplicity; such interactions could actually be used to assist the folding and will be the subject of future study.

The definition of the long-range order parameter $`S`$, Eq. (114), must be modified to reflect the structure of the anticipated collapsed state. The partial contributions to $`S`$ of the two helical segments are now combined with opposite signs, and contributions from the bridging region ignored; this provides a reasonably sensitive, but unambiguous, measure of folding success. Fig. 12 shows how both the modified $`S`$ and the energy vary with time, for $`L=78`$ chains (corresponding to a folded state consisting of a pair of 6-turn helices), using runs whose details are otherwise similar to those for $`L=90`$. Based on visual analysis of the behavior, a different definition of what constitutes successful folding is needed here, namely that the value of the modified $`S`$ must now exceed 0.95 for folding to be considered successful; using this criterion the success fraction was found to be 0.66.

Figure 13 shows an example of a successfully folded helix pair, an ordered conformation in which both secondary and tertiary structural elements are manifest; two-thirds of the runs ended in this state. The failure to fold properly was generally not due to defects in the individual helical segments, but because the two secondary components did not succeed in aligning correctly; an example of such an outcome is shown in Fig. 14. The attractive forces become much more effective once the helical segments have formed; this is due to the linear arrangement of the attraction sites enabling them to function cooperatively, an effect that may be reflected in real proteins with prominent secondary-structural features. As a result of the nature of the interactions and the relative interaction strengths, the helical segments form first, essentially unimpeded, and only then do they attempt to align. The failure to align here is a symptom of the absence of any driving force for bringing the helices together; there is no torsional preference in the bridging segment and the range of the inter-helix attraction is too short to be felt if the helical segments are well separated. Changes to either or both these aspects of the potential should alter the behavior, but care is required to avoid hindering the helix formation process in any way.

## V Conclusions

The present paper has focused on both methodology and results. A formalism developed for the dynamics of robotic manipulators and other coupled mechanical systems – that provides a convenient and direct representation of the dynamics of bodies connected by rigid links with restricted degrees of freedom – has been utilized in a polymer context, with the clear implication that existing methods based on geometric constraints may be redundant in many instances. Since the treatment also involves dealing with rigid-body dynamics, a computationally more effective method than the often-used quaternion approach is also employed.

The results of an extensive series of MD simulations demonstrate that homopolymer chains with suitable torsional interactions consistently collapse into well-formed helices; the probability of localized defects being frozen into the structure depends on the cooling rate, and it can be reduced to a very low level by cooling sufficiently slowly. In order to demonstrate that the present simplified approach is relevant for protein folding, heterogeneous chains, with interactions favoring the development of antiparallel pairs of helices, were shown to produce coexisting secondary and tertiary structural features.

The order parameters introduced to quantify the degree of folding were tailored to capture the structural order present in the final state of the polymer. To study the details of folding pathways, other order parameters (or reaction coordinates) that capture features present in the intermediate states, but not necessarily in the final state, could be defined. In the twin-helix case, for example, a simple sum of the absolute values of $`S`$, evaluated separately for each helix-forming segment of the chain, might prove useful, since this quantity reaches its maximum upon completion of secondary structure formation, and is not seriously affected by subsequent rearrangement at the tertiary level.

The apparent success of the MD approach to chain folding used here is important for another reason. The widely-cited Levinthal “paradox” Pande et al. (2000) implies that since the number of states accessible to a protein grows exponentially with residue count, the time required for even a small protein to seek out its native state is, for practical purposes, infinite. Since nature does not suffer from this problem, the implication is that substantial portions of the folding process occur along certain well-characterized pathways; thus the molecules do not really wander almost aimlessly through conformation space, and hence there is no paradox. In order to begin to simulate such processes it is necessary to resort to a computationally efficient model, with realistic dynamics and a unique but readily determined low-energy “native” state; this is precisely what has been accomplished in the present work.

The type of model introduced here provides a starting point for exploration in several directions. While the interactions were weighted to construct the secondary helix structure prior to forming features at the tertiary level, a change in the relative strength of the interactions would allow aspects of both levels of organization to appear concurrently. Chains could be designed to fold into other idealized compact structures, such as the packed cube used in some lattice studies Dill et al. (1995); Shakhnovich (1997), or sheetlike conformations that also represent important secondary structure components; furthermore, packed states with different degrees of accessibility could provide useful information on how this feature influences folding success. The interactions can be modified and new types of interactions added; polymer topology can be changed by the addition of side chains corresponding to residues with extended structure. Common to all these enhancements is that the model must always be designed with a known lowest-energy state, and in this respect the approach differs from many other types of protein simulation. While models of this kind are perhaps limited in the kinds of questions they can address, there are more than enough issues requiring attention where they can prove helpful.

In a sense, the role played by such highly simplified models is analogous to the Ising model of ferromagnetism Huang (1963); while it is not usually claimed that an Ising spin system accurately represents a real magnetic material (or, for that matter, any other kind of real physical system, when used for other kinds of problems such as lattice gases) it does however capture a great deal of the essence of the problem, to an extent that the study of Ising and related models has resulted in important advances, both for spin systems in particular, and for statistical mechanics and critical phenomena in general. Proteins can also be modeled with a high level of detail and specificity, but the tradeoff is that only short trajectory segments can be followed with an investment of a reasonable amount of computing effort; hopefully, extensive studies of simplified polymer models of the kind examined here, in which the design is tailored to reproduce certain generic aspects of macromolecular behavior, will achieve greater popularity as their usefulness becomes established.

###### Acknowledgements.

This work was partially supported by the Israel Science Foundation.


text:
# History-dependent Maximum Static Friction and Rectification of a System with a Few Particles in a Periodic Field

## Abstract

History-dependent maximum static friction is investigated using simple models with two and three particles in one-dimensional periodic potentials. In some situations, these systems possess two values of the maximum static friction, with that actually realized being determined by the direction of the slippage relative to the direction of the previous slippage. In particular, the maximum static friction for slippage in a given direction is smaller when the previous slippage was in the same direction than when it was in the opposite direction. Owing to this property, a particle can continue to move in a direction even in the case that it is subject to an external force whose direction change in time as in a ratchet system which is regarded as one of the simplest model of molecular motors. PACS number(s):

Static and dynamic frictions and the transition between them are universal phenomena commonly observed at the surfaces of macroscopic objects. There is a rich variety of such phenomena, and their characteristics depend on the physical properties of the surfaces as well as any lubricants that might exist between them. In some cases, frictions can depend on the system history . For example, the maximum static friction can depend on how long the two surfaces have been in contact . Also, in some case, dynamic friction with increasing slippage velocity is larger than that with decreasing velocity.

In this paper, the existence of history-dependent maximum static friction is investigated using simple systems. First, a system consisting of only two particles in a spatially periodic field is considered. In this system, the maximum static friction depends on the direction of the previous slippage. We show that this property can cause the system to exhibit rectification behavior. We also show that similar phenomena are observed in a system consisting of three particles in a spatially periodic field.

We study a one dimensional system containing two particles in which both particles are subject to a spatially periodic potential and one particle is subject to an external force $`F_{ex}`$ (Fig. 1). The motion of these two particles is described by the over-damped equations,

$$\dot{x_1}=c_1c_2\mathrm{sin}(2\pi (x_1x_2))+c_1c_p\mathrm{sin}(2\pi x_1)+F_{ex}$$
(1)
$$\dot{x_2}=c_1c_2\mathrm{sin}(2\pi (x_2x_1))+c_2c_p\mathrm{sin}(2\pi x_2)$$
(2)
where $`x_i`$ is the $`i`$th particle’s position. We set $`c_1=c_p=1`$ and $`c_2=c>0`$, and we consider the system to be defined for $`0x_i<1`$, with periodic boundary conditions. This system belongs to a class of coupled phase oscillators systems that have been studied extensively.

The above system can be regarded as a simplified model of a physical system in which thin lubricants is spread uniformly between two objects with bumpy surfaces, and an external force acts on one of the two objects. In this situation, the motion of the contact points of the two objects and the motion of the lubricant particles are, approximately, modeled by the motion of the first and second particle, as described by Eqs. (1) and (2), respectively (Fig. 1). The condition $`c_i>0`$ means that, as the interactions in the system, we consider only the effects of repulsive forces like the excluded volume effect which plays important roles in liquids, solids and gels.

In the following, we report the results of simulations employing above system whose purpose is to determine the relationships between the maximum static friction and observed microscopic states. We use $`F_{ex}=(1\mathrm{cos}(t/T))/2`$ with large $`T`$, so that the external force varies smoothly and very slowly. The data points ($`\times `$) in Fig. 2 indicate values of the maximum static friction $`R`$, plotted as a function of $`c`$ $`(0c1)`$ for two particle configurations, $`x_1^b>x_2^b`$ and $`x_1^b<x_2^b`$, before slippage. Here, $`R`$ is defined as the maximum value of $`|F_{ex}|`$ for which the first particle remains stuck as $`|F_{ex}|`$ increases, and $`x_i^b`$ is the $`i`$th particle’s position at the most recent time at which $`F_{ex}=0`$. This $`c`$ \- $`R`$ relation is divided into three regions: I) for $`c<\widehat{c}_{crit^1}0.59`$, $`R`$ is decreasing function of $`c`$ for all $`x_1^b`$, $`x_2^b`$, II) for $`\widehat{c}_{crit^1}<c<\widehat{c}_{crit^2}0.66`$, $`R`$ has two possible values, one realized for $`x_1^b>x_2^b`$ and one for $`x_1^b<x_2^b`$, and III) for $`c>\widehat{c}_{crit^2}`$, $`R`$ is an increasing function of $`c`$ for all $`x_1^b`$, $`x_2^b`$.

When $`|F_{ex}(t)|`$ decreases from some initial value greater than $`R`$, the first particle ceases slipping at a value of $`|F_{ex}|=R_{stop}`$. Through our simulations, we found relations between $`R_{stop}`$ and $`R`$ for each of the above described cases: $`R_{stop}=R`$ in case I), $`R_{stop}<R`$ in case III) as shown in Fig. 2, and $`R_{stop}=R_{smaller}`$ in case II), where $`R_{smaller}`$ is the smaller of the two values of $`R`$.

Figure 3 displays typical temporal evolutions of each particle’s velocity and position for cases I) and III), with $`F_{ex}`$ as given above. Here, in (a) $`c=0.2`$ and $`x_1^b<x_2^b`$, in (b) $`c=0.2`$ and $`x_1^b>x_2^b`$, in (c) $`c=0.8`$ and $`x_1^b<x_2^b`$, and in (d) $`c=0.8`$ and $`x_1^b>x_2^b`$. In this figure, the gray curves represent position and velocity of the first particle, and the black curves represent those of the second. In the situation depicted in Fig. 3 (a) and (b), $`x_1>x_2`$ always holds just before the slippage of the first particle, independently of $`x_1^b`$ and $`x_2^b`$. This is due to the fact that if $`F_{ex}`$ is small, for sufficiently small $`c_2`$, it can be the case that the first particle is not able to cross the potential barrier, while it is able to cross the second particle. In fact this is the case in the situations described by Figs. 3(a) and (b). Contrastingly, in Figs. 3 (c) and (d), $`x_1<x_2`$ always holds just before the slippage of the first particle, independently of $`x_1^b`$ and $`x_2^b`$. This follows from the fact that if $`F_{ex}`$ is small, for sufficiently large $`c_2`$, it can be the case that the first particle cannot cross the second particle, while it can cross the potential barrier. From these considerations, it is clear why $`R`$ has only a single value for each $`c`$ value in cases of I) and III).

By considering the balance equations obtained by setting $`\dot{x_1}=0`$ and $`\dot{x_2}=0`$ in Eqs. (1) and (2) and choosing the particle configuration before the slippage, we can obtain $`c`$-$`R`$ curves plotted in Fig. 2. For example, the curve of $`R`$ obtained as the maximum values of $`F_{ex}`$ with $`\dot{x_1}=0`$ and $`\dot{x_2}=0`$ under the condition $`x_1>x_2`$ is consistent with the numerical result for $`c<\widehat{c}_{crit^2}`$, where $`R=[\mathrm{sin}(2X)+c\mathrm{sin}(X)]`$ with $`X=\mathrm{arccos}[(c(c^2+32)^{1/2})/8]`$. The line $`R=c`$ is consistent with the numerical results for $`c>\widehat{c}_{crit^1}`$, which is obtained from $`R=[c\mathrm{sin}(2\pi (x_1x_2))+\mathrm{sin}(2\pi x_1)]|_{x_1=0.5,x_2=0.75}`$. Here $`x_1=0.5`$ and $`x_2=0.75`$ corresponds to the values for which the force of the second particle on the first under the condition $`x_1<x_2`$ is maximal.

In case II), the maximum static friction is determined by the history, which direction did the first particle slip previously. This can be understood as follows. Figure 4 displays typical temporal evolutions of each particle’s velocity and position for case II). (Here, $`c=0.63`$ and in (a) $`x_1^b<x_2^b`$ and in (b) $`x_1^b>x_2^b`$.) In contrast to cases I) and III), here, the particles do not cross as $`F_{ex}`$ is increased until the first particle starts to slip. We also found that slippage begins at a later time for $`x_1^b<x_2^b`$ than for $`x_1^b>x_2^b`$ . This means that $`R`$ depends on the particle configuration: $`R`$ for $`x_1^b<x_2^b`$ is larger than that for $`x_1^b>x_2^b`$. Moreover, the state with $`x_1>x_2`$ is realized after slippage of the first particle whether $`x_1^b>x_2^b`$ or $`x_1^b<x_2^b`$. (Because of the symmetry of this system, $`x_1<x_2`$ is necessary realized after slippage of the first particle if $`F_{ex}<0`$.) This means that $`R`$ for slippage in the same direction as the previous slippage is always smaller than that for the slippage in the opposite direction.

If $`F_{ex}>0`$ and it varies slowly and smoothly in time, in case II), $`x_1>x_2`$ is realized after slippage. We now explain the reason for this. Figure 4 (c) plots the temporal evolutions of each particle’s velocity and position for $`c=0.63`$ with $`F_{ex}=0.583`$, which is slightly larger than $`R_{smaller}`$. If $`F_{ex}`$ is constant and $`|F_{ex}|>R`$, the first particle slips and the second particle oscillates periodically. As shown in Fig. 4 (c), the time required to switch from $`x_1>x_2`$ to $`x_1<x_2`$ is much longer than that to switch from $`x_1<x_2`$ to $`x_1>x_2`$ when $`F_{ex}`$ is slightly larger than $`R_{smaller}`$. If $`c`$ is not large, the amplitude of the second particle’s oscillation is not large. In such a situation, $`x_2`$ cannot reach such a value ($`x_20.75`$) that the force of the second particle on the first is large enough to balance $`F_{ex}`$ at $`x_1<x_2`$. Hence, the system cannot remain in a state with $`x_1<x_2`$ for an extended time if the first particle crosses over the potential barrier. in contrast to case III) depicted in Figs. 3 (c) and (d).

Also in case III), the second particle cannot reach at $`x_20.75`$ when $`|F_{ex}|`$ decreases from a value greater than $`R`$. However, in this case, the $`x_2`$ can reach such a value that the force of the second particle on the first at $`x_1<x_2`$ balances $`F_{ex}`$ with $`|F_{ex}|=R_{stop}`$, in contrast to case II).

Since, in case II), the static maximum friction to the direction as the direction of previous slippage is smaller than that to the opposite direction, for this case, this system exhibits rectification behavior. Recently, rectification behavior has been observed, for example, in some ratchet systems that are simple models of molecular motors. In these models, with an external force that favors neither direction, particles can move in only one direction, which is completely determined by the shape of the potential. Now, to demonstrate that our model can exhibit similar behavior, we consider the case $`c=0.63`$ and $`F_{ex}=0.61\mathrm{sin}(t/T)`$, with sufficient large $`T`$. Thus, in this case, the external force acts symmetrically in both positive and negative directions. Here, we choose the value $`0.61`$ because it is halfway between the two values of $`R`$ for $`c=0.63`$. For this system, as shown in the left-hand side of Fig. 4 (d), the first particle can move only in one direction, even though $`F_{ex}`$ alternates between positive and negative values. However, if a sufficient strong force is applied instantaneously to the first particle in the direction opposite to its motion, its motion can reverse direction \[see the right-hand side of Fig. 4 (d)\]. This means we can control the slippage direction in the system. This is not possible in the ratchet systems proposed as models of molecular motors.

The phenomena discussed above can also be observed in systems with a larger number of degrees of freedom. As an example, we consider a system that consists of three particles in a spatially periodic field, with the motions of each particle obeying the equation,

$$\dot{x_i}=\underset{ij}{}c_ic_j\frac{1+\mathrm{cos}(2\pi (x_ix_j))}{2}\mathrm{sin}(2\pi (x_ix_j))+c_ic_p\frac{1+\mathrm{cos}(2\pi x_i)}{2}\mathrm{sin}(2\pi x_i)+\delta _{i,1}F_{ex}.$$
(3)
The characteristic length of the interactions between particles in this system is shorter than that of the system described by (0.1) and (0.2). Again, we consider the system to be defined in the region $`0x_i<1`$ and use periodic boundary conditions. In the case that $`c_1=c_2=1`$ and $`c_3=c_p=c`$, behavior similar to that exhibited in cases II) for the two-particle system is observed over a wide range of values of $`c`$ with $`F_{ex}=F(1\mathrm{cos}(t/T))/2`$, for $`F>0`$ and sufficiently large $`T`$.

The solid curves in Fig. 5 represent $`R`$ as a function of $`c`$ $`(0.1c1)`$. As shown, $`R`$ takes two values for each value of $`c`$, depending on the relationships among the $`x_i^b`$. Figures 6 (a), (b) and (c) displays typical temporal evolutions of each particle’s velocity and position with $`c=0.5`$ for two particle configurations before slippage; the (1,2,3) configuration in which the particles are arranged in the order $`1`$, $`2`$, $`3`$ with respect to the direction of $`F_{ex}`$ (the direction of increasing $`x`$ in this case), as in (a), and the (1,3,2) configuration, as in (b) and (c). Here, $`F=0.5`$ in (a) and (b), and $`F=0.36`$ in (c). By comparing (a) and (b), we find that $`R`$ for the (1,3,2) configuration is larger than that for the (1,2,3) configuration. Also it is seen that the (1,2,3) configuration is always realized after slippage in both cases considered in (a) and (b). We are thus led to the conclusion that, as in the system with two particles, under a certain condition, $`R`$ for the direction of the previous slippage is always smaller than $`R`$ for the direction opposite to the previous slippage. For the three particle system, this condition is that $`F`$ be larger than a particular value, which we discuss below.

As stated above, the (1,2,3) configuration is apparently always realized after slippage in the situations considered in Figs. 5 (a) and (b). If $`F`$ is not too large, however, the (1,3,2) configuration can be preserved upon slippage. In fact, this is the case for the situation considered in Fig. 6 (c). The dotted curve in Fig. 5 represent critical values of $`F_{ex}`$, below which the particle configuration is preserved upon slippage. Thus, in this system, the maximum strength of the external force at the slippage determines whether or not there is a direction dependence of $`R`$. We found that with a properly chosen $`F_{ex}`$, this system too can exhibit rectification behavior.

In this paper, we have investigated the history dependence of the maximum static friction using simple systems consisting of two and three particles in a one-dimensional periodic potential. In these system, we found that, in some cases, the maximum static friction can depend on the direction of the slippage, being smaller than this direction is the same as that of the previous slippage. By this property, a particle in this system can continue to slip along a single direction even in the case that the direction of the external force change in time. This behavior is similar to that seen in ratchet systems regarded as one of simplest models of molecular motors. Analytical study of the systems considered here and further investigation of systems with three or more particles will be carried out in the future. The construction of a model that can realize more types of memory effects, for example aging of friction, is also an important future problem.

The author is grateful to K. Kaneko, K. Fujimoto, K. Sekimoto and H. Matsukawa for useful discussions. This research was supported in part by a Grant-in-Aid for JSPS Fellows (10376).