New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Nonlinear relaxation dynamics in elastic networks and design principles of molecular machines

Communicated by Gerhard Ertl, Max Planck Society for the Advancement of Science, Berlin, Germany, April 4, 2007 (received for review January 23, 2007)
Abstract
Analyzing nonlinear conformational relaxation dynamics in elastic networks corresponding to two classical motor proteins, we find that they respond by well defined internal mechanical motions to various initial deformations and that these motions are robust against external perturbations. We show that this behavior is not characteristic for random elastic networks. However, special network architectures with such properties can be designed by evolutionary optimization methods. Using them, an example of an artificial elastic network, operating as a cyclic machine powered by ligand binding, is constructed.
Understanding design principles of singlemolecule machines is a major challenge. Experimental and theoretical studies of proteins, acting as motors (1–5), ion pumps (6–8), or channels (6, 9), and enzymes (10–14), show that their operation involves functional conformational motions (see ref. 15). Such motions are slow and cannot therefore be reproduced by full molecular dynamics simulations. Within the last decade, approximate descriptions based on elastic network models of proteins have been developed (16–21). In this approach, structural elements of a protein are viewed as identical point particles, with two particles connected by an elastic string if the respective elements lie close enough in the native state of the considered protein. Thus, a network of elastic connections corresponding to a protein is constructed. So far, the attention has been focused on linear dynamics of elastic networks, characterized in terms of their normal vibrational modes. It has been found that ligandinduced conformational changes in many proteins agree with the patterns of atomic displacements in their slowest vibrational modes (refs. 22–25 and see also refs. 26 and 27), even though nonlinear elastic effects must become important for large deviations from the equilibrium (28, 29). The focus of this article is on nonlinear relaxation phenomena in elastic networks seen as complex dynamical systems.
Generally, a machine is a mechanical device that performs ordered internal motions that are robust against external perturbations. In machines representing single molecules, energy is typically supplied in discrete portions, through individual reaction events. Therefore, their cycles consist of the processes of conformational relaxation that follow after energetic excitations. For a robust machine operation, special nonlinear relaxation dynamics is required. We expect that, starting from a broad range of initial deformations, such dynamical systems would return to the same final equilibrium state. Moreover, the relaxation would proceed along a well defined trajectory (or a lowdimensional manifold), rapidly approached starting from different initial states and robust against external perturbations. These attractive relaxation trajectories would define internal mechanical motions of the machine inside its operation cycle.
This special conformational relaxation dynamics has been confirmed in our study of the elastic networks of two protein motors (F_{1}ATPase and myosin). On the other hand, our control investigation of random elastic networks has shown that relaxation patterns in random elastic networks are typically complex and qualitatively different from those of protein motors. Actual proteins with specific architectures allowing robust machine operation may have developed through a natural biological evolution, with the selection favoring such special dynamical properties. In a model study, we have demonstrated that artificial elastic network architectures possessing machinelike properties can be designed by running an evolutionary computer optimization process. Finally, an example of an artificially designed elastic network that operates like a machine powered by ligand binding has been constructed.
Elastic Network Models
The considered elastic networks consist of a set of N identical material particles (nodes) connected by identical elastic strings (links). A network is specified by indicating equilibrium positions of all particles. Two particles are connected by a string if the equilibrium distance between them is sufficiently small. The elastic forces, acting on the particles, obey Hooke's law and depend only on the change in the distances between them. In the overdamped limit (18), the velocity of a particle is proportional to the sum of elastic forces applied to it. If R _{i} ^{(0)} are equilibrium positions of the particles and R _{i} (t) are their actual coordinates, the dynamics is described by: where A is the adjacency matrix, with the elements A_{ij} = 1, if R _{i} ^{(0)} − R _{j} ^{(0)} < l _{0}, and A_{ij} = 0 otherwise. The dependence on the stiffness constant of the strings and the viscous friction coefficient of the particles is removed here by an appropriate rescaling of time.
The dynamics of elastic networks is nonlinear, because distances R
_{i}
− R
_{j}
 are nonlinear functions of the coordinates R
_{i}
and R
_{j}
. Close to the equilibrium, equations of motion can however be linearized, yielding:
for small deviations r
_{i}
= R
_{i}
− R
_{i}
^{(0)}. These equations can be written as ṙ
_{i}
= − Σ
_{j}
Λ
_{ij}
r
_{j}
, where Λ is a 3N × 3N linearization matrix. In the linear approximation, relaxation motion is described by a sum of independent exponentially decaying normal modes:
with λ_{α} and e
_{i}
^{(α)} representing nonzero eigenvalues and the respective eigenvectors of the matrix Λ. It should be noted that the same eigenvalues determine vibration frequencies of the network, ω_{α} ∼
Results
Two Motor Proteins.
As an example, Fig. 1 shows the elastic network of the single βsubunit of the molecular motor F_{1}ATPase (Protein Data Bank ID code 1H8H, chain E) (30). Each node corresponds to a residue in this protein (the total number of nodes is N = 466). In Fig. 1 b, a set of conformational relaxation trajectories, obtained by numerical integration of the nonlinear elastic model (see Methods and Eq. 1 ) of this macromolecule, is displayed. To track conformational motions, three network nodes (1, 2, and 3) were chosen, and pair distances u _{12}, u _{13}, and u _{23} were determined during the relaxation process. Thus, each conformational motion was represented by a trajectory in a 3D space, where coordinates were normalized deviations Δu_{ij} /u _{ij} ^{(0)} from the equilibrium pair distances u _{ij} ^{(0)} = R _{i} ^{(0)} − R _{j} ^{(0)} (u_{ij} = R _{i} − R _{j}  and Δu_{ij} = u_{ij} − u _{ij} ^{(0)}). Each trajectory begins from a different initial conformation obtained by applying random static forces to all network nodes (see Methods). Trajectories starting from various initial conditions soon converge to a well defined relaxation path leading to the equilibrium state. This path corresponds to a slow motion of the network group, including label 1, with respect to the rest of the molecule. The protein is soft along such a path: by applying static forces of the same magnitude, one can stretch it by 30% along the path, as compared with the length changes of only a few percent when the forces were applied in other directions.
Relaxation behavior in the nonlinear elastic network (N = 793) of another classical molecular motor, myosin (single heavy chain; Protein Data Bank ID code 1KK8, chain A) (31), is displayed in Fig. 2. In contrast to the βsubunit of F_{1}ATPase, this elastic network possesses an attractive 2D manifold (a plane). The network is extremely stiff for deformations in the directions orthogonal to this plane. By applying static forces of the same magnitude, one can only induce relative deformations of ≈10^{−3} along such orthogonal directions, as compared with the relative deformations of the order of 10^{−1} for the directions within the plane. To characterize the temporal course of relaxation, one relaxation trajectory (displayed in green in Fig. 2 b and c) and subsequent positions at different time moments along the trajectory are indicated in Fig. 2 d. The trajectory rapidly reaches the plane and then the relaxation motion becomes much slower (with the characteristic time scales larger by a factor 10–100). A similar behavior is characteristic for other relaxation trajectories, starting from different initial conditions. All recorded trajectories returned to the equilibrium state, and no metastable states were encountered starting from the chosen initial conditions.
Nonlinear effects were essential in the relaxation dynamics starting from large arbitrary initial deformations considered here. Remarkably, the observed relaxation patterns are nonetheless qualitatively in agreement with the normal mode analysis. Both motors possess a group of soft modes separated by a gap from the rest of the spectrum (eigenvalue spectra of elastic networks of these proteins are shown in Fig. 3). The two soft modes of myosin define the attractive plane seen in the relaxation pattern of its elastic network (Fig. 2). The elastic network of F_{1}ATPase has three soft modes. However, one of the soft modes has a relaxation rate that is smaller than the other two modes. Therefore, the pattern of relaxation trajectories looks here like a thick 1D bundle. Specific ligandinduced conformational changes in F_{1}ATPase and myosin were previously shown to have strong overlaps with patterns of deformation in slow vibrational modes (24).
Random and Designed Elastic Networks.
A control study of nonlinear relaxation phenomena in random elastic networks was performed. Such networks were obtained by taking a relatively short chain of N = 64 and randomly folding it in absence of energetic interactions (see Methods). After that, all particles separated by short enough distances were connected by identical elastic links. Fig. 4 shows relaxation patterns for two typical random elastic networks (using the same procedure for generation of initial conditions and for tracking of conformational relaxation as in Figs. 1 and 2). Relaxation patterns in random networks are clearly different from those of the considered motor proteins. Random networks possess many (meta)stable stationary states, with relaxation trajectories often ending in one of them instead of going back to the equilibrium conformation. The linear normal mode description holds in such networks only in close proximity of the equilibrium state.
Thus, elastic networks of motor proteins are special. Their equilibrium conformation has a big attraction basin. Starting from an arbitrary initial deformation, relaxation dynamics is soon reduced to a lowdimensional attractive manifold. Within its large part, the dynamics is approximately linear and determined by a few soft modes. Proteins with such special dynamical properties, essential for their functions, may have emerged through a biological evolution. Below, we show that artificial elastic networks with similar properties can be constructed by running an evolutionary optimization process based on a variant of the Metropolis algorithm.
The evolutionary optimization technique is described in Methods. For each network, spectral gap g = log_{10}(λ_{2}/λ_{1}) is defined as the logarithm of the ratio between the relaxation rates λ_{2} and λ_{1} of its two slowest normal modes. If a substantial gap is present, the slowest mode has the relaxation rate that is much smaller than that of the other modes; therefore, the longtime relaxation in the linear regime would be dominated by this soft normal mode. The used evolutionary optimization process maximizes the spectral gap g of the evolving networks. Beginning with an initial random network, we applied structural mutations and determined the difference Δg = g′ − g of the gaps before and after a mutation. If the gap was increased (Δg > 0), the mutation was always accepted. If Δg < 0, the mutation was accepted with the probability P = exp(Δg/θ), where θ is the effective optimization temperature. This procedure was applied iteratively, generating an evolution started from an initial network. ^{†}
Networks with soft normal modes and large spectral gaps were thus constructed. A typical network with a large gap and its relaxation pattern are shown in Fig. 5 a and b [for other examples, see supporting information (SI) Fig. 7)]. The presence of a large gap has a strong effect on the nonlinear relaxation properties in such systems. They possess well defined long paths with slow conformational motion leading to the equilibrium state. There are only a few (meta)stable states and these states usually lie on an attractive relaxation path, so that a small barrier is encountered when moving along the path. The opposite behavior with complex relaxation patterns and a high number of (meta)stable conformations was found for a set of “failed” networks where gaps were small and could not be significantly increased through the evolution (see SI Fig. 8). Our analysis shows that the designed networks can be viewed as consisting of rigid blocks connected by soft joints; they are able to perform some large conformational changes accompanied by only small local deformations (see SI Fig. 9).
Spectral gaps of the networks, designed by using this optimization method, are much larger than those characteristic for motor proteins (compare Fig. 3). To improve the agreement, additional “reverse” evolution was subsequently applied to the designed networks, with the selection pressure aimed to decrease the gap (see Methods). While the gap was rapidly reduced, the global relaxation pattern was changing only much more slowly with structural mutations and retained characteristic features of the networks with large gaps. Fig. 5 c displays the network, obtained by applying such reverse evolution (with only five subsequent mutations) to the original network shown in Fig. 5 a. Although the spectral gap has been reduced from 9.53 to 1.22, the principal structure of the network is retained, with the mutations mostly affecting only the hinge region. The relaxation pattern of the constructed network (Fig. 5 d) reveals an attractive path leading to the equilibrium state (with another stable state lying on it). Remarkably, the linear normal mode description of relaxation dynamics holds in such networks within a much larger domain around the equilibrium state. We conclude that, by running an evolutionary optimization process, artificial elastic networks approaching conformational relaxation properties of real motor proteins can be constructed.
An Artificial Machine Network.
Using such designed networks, we proceed to construct nonlinear elastic systems that can be viewed as prototypes of a machine powered by ligand binding. The network shown in Fig. 6 performs cyclic hinge motions, caused by binding and detachment of a ligand. To obtain it, an initial network with two distinct parts (dense clusters), which were only loosely connected, was prepared. This initial network was characterized by a small spectral gap. By running evolutionary optimization, the architecture of the network was modified, so that it developed a large spectral gap while maintaining its special structure. The equilibrium conformation of the finally obtained network is shown in Fig. 6 Lower Right. The cycle began with binding of a ligand (snapshot A, t = 0). The ligand was an additional particle that could form elastic connections with the existing network nodes. Binding of a ligand was modeled by placing the particle at a location chosen in the hinge region and creating elastic links with the three near nodes (for details, see Methods). When such new links were introduced, they were in a deformed (i.e., stretched) state and a certain amount of energy was thus supplied to the system. Therefore, the networkligand complex was not initially in the equilibrium state and conformational relaxation to the equilibrium had to occur. Within a short time, deformations spread over a group of near elastic links (B, t = 200), producing a small change in the network configuration. After that, a slow hinge motion of the networkligand complex took place (C, t = 4,000) until the equilibrium state of the complex was reached (D, t = 20,000). In this state, the ligand was removed. At that moment, the elastic network was in a configuration far from equilibrium and its relaxation set on. Within a short time, the network reached the path (E, t = 20,200) along which its slow ordered relaxation proceeded (F, t = 30,000), until the equilibrium conformation (A) was reached. Then, a new ligand could bind at the network, and the cycle was repeated. For visualization of conformational motions inside the network cycle, see SI Movie 1.
To monitor conformational motions, three labels were attached to the network, and distances u _{12} and u _{13} between them were recorded. Fig. 6 Lower Left shows trajectories in the plane (Δu _{12}, Δu _{13}) that correspond to the operation of this prototype machine in the presence of noise and with stochastic binding of ligands. As seen in Fig. 6, the trajectories of the forward and back hinge motions are different, but both of them are well defined. The applied noise (see Methods) is only weakly perturbing them, which means that both the free network and the network–ligand complex have narrow valleys with steep walls, leading to their respective equilibrium states. Binding of a ligand is possible in a relatively broad interval of conformations near the equilibrium and therefore there is a dispersion in the transitions from the upper to the lower branch of the cycle. The operation of this machine is essentially nonlinear and unharmonic effects are strong. Characteristic rates of slow relaxation motions in this network differ by a factor between 10 and 1,000 from the rate of the fast conformational transition after the ligand binding.
Discussion
In conclusion, we have shown that motor proteins possess unique dynamical properties, intrinsically related to their functioning as machines. We have also demonstrated that artificial elastic networks with similar properties can be constructed by evolutionary optimization methods. To verify these theoretical predictions, special singlemolecule experiments monitoring conformational relaxation after arbitrary initial deformations can be performed. An example of an elastic machinelike network powered by ligand binding is presented. Using such designed networks, fundamentals of molecular machine operation, including energetic aspects, the role of thermal noise and hydrodynamic interactions, can be discussed (32). Comparing the behavior of such designed networks with that of real molecular machines, better understanding of what is general and what is specific for a particular protein can be gained. Moreover, our analysis provides a systematic approach for the design of proteins with prescribed (programmed) robust conformational motions. Not only proteins, but also atomic clusters can be described by elastic network models (33). Similar methods can further be used for engineering of nonprotein machinelike nanodevices (see ref. 34).
Methods
The elastic networks for protein structures in Figs. 1 and 2 are constructed from the structural data of F_{1}ATPase in the ATPanalog state (Protein Data Bank ID code 1H8H) chain E (30) and myosin in the ATPanalog state (Protein Data Bank ID code 1KK8) chain A (31) with a cutoff distance l _{0} = 10 Å. We place a node at the position of each αcarbon atom and connect nodes lying within the cutoff distance by a link. Stiffness constants of all links and friction coefficients for all particles are equal.
To generate random networks, a chain of N nodes is taken and folded randomly in the 3D space. Each next node is positioned at random, in such a way that its distance l from the preceding node satisfies the condition l _{min} ≤ l ≤ l _{max}. Moreover, the next node should be separated by at least the distance l _{min} from all previous nodes. Examining the folded chain, all pairs of nodes (i, j) with distances <l _{0} are noticed (l _{0} > l _{max}) and connected by elastic strings (hence, not only all neighbors in the string get connected, but also those nodes that come by chance close one to another in the folded conformation).
The iterative optimization process consists of a sequence of structural mutations followed by selection. To perform a mutation, a node is taken at random and its equilibrium position is changed. The new equilibrium position is chosen with equal probability within a sphere of radius d around the old equilibrium position. To preserve the backbone chain, we additionally require that, after a mutation, the distances of the node from its both neighbors in the chain lie within an interval from l _{min} to l _{max}; other pair distances should not be shorter than l _{min} either. The new graph of connections after a mutation is constructed by examining the pair distances between all nodes and linking the nodes separated at equilibrium by a distance shorter than l _{0}.
Sometimes, networks allowing internal free rotations (and, thus, additional zero eigenvalues) may be obtained. To distinguish zero eigenvalues numerically, we adopt a numerical cutoff δλ = 10^{−12} and assume the eigenvalues < δλ to be zero. For the networks without internal rotation modes, that is, if the number N_{nz} of nonzero modes is equal to 3N − 6, we proceed in each iteration step as described. When such modes are present, a mutation is always accepted when it decreases the number of zero eigenvalues (N′_{zm} < N_{zm} ) and accepted with the probability P = exp[−(N′_{zm} − N_{zm)} /θ] otherwise. In this study, we consider networks with N = 64 nodes. The parameter values d = l _{min} = 3.4, l _{max} = 4.2, l _{0} = 8.0, and θ = 0.1 are used.
Most of the networks after optimization (97.3%) had a large gap g > 3, whereas such gaps were rare in the initial networks (1.9%). Here, only networks without internal rotations were counted (1,511 random networks of 30,000 trials and 2,346 selected networks of 2,500 trials).
In the reverse evolution used to generate special networks with smaller gaps (such as shown in Fig. 5 c), the same evolution algorithm is used with the replacement of g by −g; the effective temperature is θ = 0.01 in these simulations.
Around the equilibrium conformation, relative changes p _{ij} ^{(α)} in the distances u_{ij} between nodes i and j in a relaxation mode α are calculated as p _{ij} ^{(α)} = ∂u_{ij} /∂X_{α} = e _{ij} ^{(α)}·u _{ij} ^{(0)}/u _{ij} ^{(0)}, where e _{ij} ^{(α)} = e _{i} ^{(α)} − e _{j} ^{(α)}, u _{ij} ^{(0)} = R _{i} ^{(0)} − R _{j} ^{(0)} and u _{ij} ^{(0)} = u _{ij} ^{(0)} (see also Eq. 3 ). Here, the eigenvectors are normalized as Σ _{i} e _{i} ^{(α)}^{2} = 1. In SI Fig. 9, statistical distributions of relative changes p _{ij} ^{(1)} in the slowest relaxation mode (α = 1) in ensembles of elastic networks are shown.
For trajectory visualizations, three labels (1, 2, and 3) are attached to a network and distances u _{12}, u _{13}, and u _{23} between them are recorded. The labels are chosen in such a way that p _{12} ^{(1)}, the relative distance change in the slowest relaxation mode, is maximal between the nodes labeled as 1 and 2; then the third node is chosen for which p _{13} ^{(2)}, in the second slowest relaxation mode, is maximal between the nodes labeled as 1 and 3 (there are two choices and we have selected one of them).
To prepare initial conditions when relaxation patterns are considered, we apply randomly chosen static forces F _{i,static} and wait for a certain time τ _{s} . These random forces are chosen in such a way that their total magnitude F_{s} = [Σ_{i=1} ^{N}F _{i,static} ^{2}]^{1/2} is fixed. Then, we remove the static forces and record the trajectory. For each network, 100 trajectories from different initial conditions are shown. In Fig. 1, F_{s} = 10, τ _{s} = 3 × 10^{4}; in Fig. 2, F_{s} = 0.1, τ _{s} = 10^{5}; and in Figs. 4 and 5 and SI Figs. 7 and 8, F_{s} = 0.1 and τ _{s} = 10^{6}.
In Fig. 6 and SI Movie 1, we initially place a ligand in the center of mass of the three ligandbinding nodes and introduce additional links to these nodes with the natural length of 1.7. At the initial moment, the links are longer than their natural lengths (i.e., they are stretched) and thus attractive forces between the ligand and the binding nodes are acting. To include fluctuations, random timedependent forces f _{i} (t) of intensity σ have been added to the equations of motion of all particles, i.e., Ṙ _{i} = F _{i,elastic} + f _{i} (t) with 〈f _{i} (t)〉 = 0 and 〈f _{i} (t)·f _{j} (t′)〉 = 6σδ _{ij} δ(t − t′). We have σ = 10^{−3} (Fig. 6 Upper and SI Movie 1) and 3 × 10^{−3} (Fig. 6 Lower Left). At t = 20,000 (which is long enough for the network to relax to the steady state with a ligand, snapshot D), we remove the ligand and the additional links. In the ligandfree state, we have assumed that the network binds with a ligand again at a constant rate (probability per unit time) ν even in nonequilibrium network configurations, provided that all distances u_{ij} between the three ligandbinding nodes satisfy the conditions Δu_{ij}  = u_{ij} − u _{ij} ^{(0)} ≤ ε _{l} , i.e., the conformation of the ligandbinding site is close enough to the initial one. Then, a ligand is placed again in the center of mass of the ligandbinding nodes and the next cycle starts. The parameter values ν = 10^{−6} and ε _{l} = 0.01 are used. Coordinates of the nodes in the constructed example of the machinelike network and positions of the three binding nodes are given in SI Data Set.
Acknowledgments
Y.T. was supported by the Japan Society for the Promotion of Science through a fellowship for research abroad (H17).
Footnotes
 To whom correspondence should be sent at the present address: Nanobiology Laboratories, Graduate School of Frontier Biosciences, Osaka University, 13 Yamadaoka, Suita, Osaka 5650871, Japan. Email: togashi{at}phys1.med.osakau.ac.jp

Author contributions: Y.T. and A.S.M. designed research; Y.T. and A.S.M. performed research; Y.T. analyzed data; and Y.T. and A.S.M. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/cgi/content/full/0702950104/DC1.

↵ † The described evolutionary optimization algorithm allows us to construct only the networks with a single soft mode. It can be, however, easily modified for the design of networks with two or more soft modes, separated by a gap from the rest of the spectrum.

Freely available online through the PNAS open access option.
 © 2007 by The National Academy of Sciences of the USA
References
 ↵
 ↵
 ↵

↵
 Kitamura K ,
 Tokunaga M ,
 Esaki S ,
 Iwane AH ,
 Yanagida T
 ↵

↵
 Gouaux E ,
 MacKinnon R
 ↵
 ↵

↵
 Perozo E ,
 Cortes DM ,
 Cuello LG

↵
 Blumenfeld LA ,
 Tikhonov AN

↵
 Adams MJ ,
 Buehner M ,
 Chandrasekhar K ,
 Ford GC ,
 Hackert ML ,
 Liljas A ,
 Rossmann MG ,
 Smiley IE ,
 Allison WS ,
 Everse J ,
 et al.

↵
 Reed LJ ,
 Hackert ML

↵
 Lerch HP ,
 Mikhailov AS ,
 Hess B

↵
 Lerch HP ,
 Rigler R ,
 Mikhailov AS
 ↵
 ↵
 ↵
 ↵

↵
 Tama F ,
 Sanejouand YH
 ↵
 ↵

↵
 Brooks B ,
 Karplus M
 ↵

↵
 Zheng W ,
 Doniach S
 ↵

↵
 Cui Q ,
 Bahar I
 ↵
 ↵

↵
 Miyashita O ,
 Onuchic JN ,
 Wolynes PG
 ↵

↵
 Himmel DM ,
 Gourinath S ,
 Reshetnikova L ,
 Shen Y ,
 SzentGyörgyi AG ,
 Cohen C
 ↵
 ↵
 ↵