pubs.acs.org/Langmuir © 2009 American Chemical Society
Contact and Impact in the Multibody Dynamics of Motor Protein Locomotion Alan P. Bowling,*,† Andre F. Palmer,‡ and Lauren Wilhelm† †
Department of Mechanical and Aerospace Engineering, The University of Texas at Arlington, Arlington, Texas 76019, and ‡William G. Lowrie Department of Chemical and Biomolecular Engineering, The Ohio State University, Columbus, Ohio, 43210 Received May 21, 2009. Revised Manuscript Received July 30, 2009
This paper presents the development of a rigid multibody dynamics approach to modeling, simulation, and analysis of motor proteins. A key element of this new model is that it retains the mass properties, in contrast to many commonly used models that do not. The mass properties are usually omitted because their inclusion yields a model with multiple time scales whose simulation requires a significant amount of time. However, the proposed model can be numerically integrated in a reasonable amount of time. Thus this approach represents a new method for treating multiple scale models. In addition, retaining the mass properties allows a detailed study of contact and impact between the protein and substrate, which is critical for protein processivity. The new model also provides insights into the characteristics of the protein and its environment, specifically, the effective damping experienced by the protein moving through its fluid environment may be quite small yielding under or critically damped motion. This conclusion runs contrary to the widely accepted notion that the protein’s motion is strictly over damped. Herein, the differences between the motion predicted by the old and new modeling approaches are compared using a simplified model of Myosin V.
1. Introduction This work examines the utility of the small mass assumption in simplifying a rigid multibody dynamics model of processive motor proteins such as Kinesin and Myosin V. They “walk” along the cytoskeleton, similar to a biped, transporting organelles and nutrient packets around the cell as needed. The goal is to develop a model for studying contact/impact between the protein and cytoskeleton. The mechanics of locomotion are studied excluding chemical kinetics.1 The contributions in this paper are as follows: • Discovery of the possibility that the motor protein may experience little effective damping as it locomotes through its fluid environment; its motion may be under damped even though the fluid’s coefficient of viscous damping is large. • An analytical approach toward simulating a multiple time scale dynamic model in a reasonable amount of time. • A motor protein model which retains the mass terms allowing a detailed study of contact and impact with respect to docking; the model gives a more realistic prediction of the protein’s locomotion. The most common modeling approach represents the protein using a small number of particles, oftentimes 2 or 4.2-4 This approach has two key problems: (1) unmodeled interparticle interactions and (2) proportionality issues between model terms. Different particles contact the cytoskeleton over time, so it is not clear which ones to model. This can be addressed by modeling all of them, as in molecular dynamic simulations. These *To whom correspondence should be addressed E-mail:
[email protected]. (1) Bolterauer, H.; Tuszynski, J. A.; Unger, E. Cell Biochem. Biophys. 2005, 42, 95–119. (2) Stratopoulos, G. N.; Dialynas, T. E.; Tsironis, G. P. Phys. Lett. A 1999, 252, 151–156. (3) Mateos, J. L. Fluctuation Noise Lett. 2004, 4, L161–L170. (4) Ciudad, A.; Sancho, J. M.; Tsironis, G. P. J. Biol. Phys. 2006, 32, 455–463.
12974 DOI: 10.1021/la901812k
large models can be simplified by grouping particles into rigid bodies.5-10 Here, this is reasonable because the protein is rigid enough to act as a lever.11,12 The resulting rigid multibody model involves multiple time scales because of the disproportional size of the small molecular masses, 0.48 ag is the total mass for Myosin V, relative to large viscous friction forces, 108 ag/ms is the viscous damping coefficient for water. Approaches for addressing this can be roughly categorized as algorithmic and analytic. Algorithmic methods involve clever restructurings of the numerical integration scheme5,9,13-15 or computational architecture16 to obtain a reasonable simulation run time. Analytic methods isolate the dynamics of the larger time scales. These include modal analysis retaining only low frequency modes,7 the method of multiple scales (MMS) retaining only low order perturbations,17-19 and the small mass assumption, which omits the mass terms.3,4,20 The relationship between the MMS and the small mass assumption was discussed in a previous
(5) Chun, H. M.; Padilla, C. E.; Chin, D. N.; Watanabe, M.; Karlov, V. I.; Alper, H. E.; Soosaar, K.; Blair, K. B.; Becker, O. M.; Caves, L. S. D.; Nagle, R.; Haney, D. N.; Farmer, B. L. J. Comput. Chem. 2000, 21, 159–184. (6) Schuyler, A. D.; Chirikjian, G. S. J. Mol. Graphics Modell. 2004, 22, 183–193. (7) Schuyler, A. D.; Chirikjian, G. S. J. Mol. Graphics Modell 2005, 24, 46–58. (8) Schwieters, C. D.; Clore, G. M. Biochemistry 2007, 46, 1152–1166. (9) Mukherjee, R. M.; Crozier, P. S.; Plimpton, S. J.; Anderson, K. S. Int. J. Non-Linear Mech. 2008, 43, 1040–1055. (10) Bowling, A.; Palmer, A. F. J. Biomech. 2009, 42, 1218–1223. (11) Purcell, T.; Morris, C.; Spudich, J. A.; Sweeney, H. L. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 14159–14164. (12) Shiroguchi, K.; Kinosita, K. Science 2007, 316, 1208–1212. (13) Jain, A.; Vaidehi, N.; Rodriguez, G. J. Chem. Phys. 1993, 106, 258–268. (14) Watanabe, M.; Karplus, M. J. Chem. Phys. 1993, 99, 8063–8074. (15) Medyanik, S. N.; Liu, W. K. Comput. Mech. 2008, 42, 569–577. (16) Alam, S. R.; Agarwal, P. K.; Hampton, S. S.; Ong, H.; Vetter, J. S. Parallel Comput. 2008, 34, 640–651. (17) Nayfeh, A. H. Perturbation Methods; John Wiley & Sons: 1973. (18) Pavliotis, G. A.; Stuart, A. M. Phys. D 2005, 2004, 161–187. (19) Coe, J. D.; Levine, B. G.; Martinez, T. J. J. Phys. Chem. A 2007, 111, 11302– 11310. (20) Chen, J. C.; Kim, A. S. Adv. Colloid Interface Sci. 2004, 112, 159–173.
Published on Web 09/09/2009
Langmuir 2009, 25(22), 12974–12981
Bowling et al.
Article
This work sufficiently illustrates the differences in behavior predicted by both models without unnecessarily complicating them. Realistic model parameters are used and there is an effort to match the experimentally observed pattern and speed of locomotion. In the next sections, the mechanical model is presented first, followed by discussions of the first- and secondorder models.
2. Mechanical Model A planar model of Myosin V, pictured in Figure 1a, is sufficient to illustrate the differences between the first- and second-order models. The features of this simplified model21 are shown in Figure 1b, which is not drawn to scale. The dynamic model can be developed using Lagrange’s equation, Kane’s method, D’Alemberts principle, or any other techniques. The mechanical model is comprised of pin-connected rigid bodies in Figure 1. The load and tail are modeled as a single rigid body without much detail because this would only complicate the model. The connections between the tail, necks, and heads are modeled as pin joints at points R, S, and P in Figure 1. Experimental studies report that the connections at point R can be modeled as a pin joint about which the necks and tail are completely free to rotate.22,12 However, the load/tail is prohibited from rotating about point R. The necks are modeled as single rigid bodies, A and B, because they are considered to be relatively stiff.21,12 Stiffness between the heads and necks is omitted. The heads’ shape is approximated by ellipses, bodies C and D, because they simplify determination of the contact forces. The binding sites on each head and the actin filament are indicated in Figure 1b. The mass centers of each body are represented by the small half-filled circles. The vectors N1 and N2 in Figure 1b define the inertial reference frame. All other reference frames are attached to the different bodies. The multibody mechanical model has the form :: _ qÞ ¼ ΣΓðq, _ qÞ AðqÞq þ bðq, ð1Þ _ qÞ ¼ bðq,
dAðqÞ D _ q_ - ðq_ T AðqÞqÞ dt Dq
ð2Þ
where q = [q1 3 3 3 q6]T contains the generalized coordinates in · Figure 1, and q_ and q. are its time derivatives of generalized velocity and acceleration. The term A(q) ∈ R66 is the mass matrix. The forces on the left of eq (1) are referred to as generalized inertia forces since they depend on mass. The forces on the right of eq (1) are referred to as generalized active forces defined as Figure 1. Picture and schematic representation of Myosin V. The schematic shows the different rigid bodies in the model. Note the difference between the point B1, defining an actin binding site, and the vector B1. PB1R is the position vector from point B1 to point R.
paper.10 The relationship between the small mass assumption and the proposed analytic approach is discussed herein. The small mass assumption yields a first-order model that cannot satisfy the initial velocity constraints which are critical to the treatment of contact and impact. The proposed second order model retains the mass properties and can satisfy both initial position and velocity conditions. In addition, the model supports the intriguing possibility that the effective viscous friction forces are quite small and may not over damp the protein’s locomotion. Langmuir 2009, 25(22), 12974–12981
ΣΓ ¼ Γcontact þ Γfriction þ Γcharge þ Γconform þ ΓBrown
ð3Þ
Γfriction ¼ -βDðqÞq_
ð4Þ
is a where β is the viscous damping coefficient and D ∈ R function of q, which transforms viscous friction forces and moments applied at the mass center of each body into generalized active forces. The terms Γcharge, Γcontact, Γconform, and ΓBrown contain forces related to Coulomb point charges, contact and impact, conformational changes, and Brownian motion. The forces acting on this model are briefly discussed in the next 66
(21) Pratt, C.; Cornely, K. Essential Biochemistry; John Wiley & Sons: 2004. (22) Dunn, A. R.; Spudich, J. A. Nat. Struct. Mol. Biol. 2007, 14, 246–248.
DOI: 10.1021/la901812k
12975
Article
Figure 2. Viscous friction forces acting on myosin V.
sections; a more detailed discussion is given in Supporting Information A. 2.1. Inertial Forces. The mass of each body is used to determine its moments of inertia. The neck regions, bodies A and B, are considered to be slender rods and the standard formulas for moments of inertia are used.23 The heads are modeled as ellipses and their moments of inertia can be found using standard formulas for an ellipsoid. The load is prohibited from rotating so its moment of inertia is not needed. The unit of mass, the attogram (ag) is chosen so that the mass values are on the order 100, and the length and time units, the nanometer (nm) and millisecond (ms), are chosen for similar reasons. These masses and inertias are contained in the mass matrix A in eq (1) which is symmetric positive definite and nondiagonal. 2.2. Viscous Friction. Each body has a force and moment applied at and about its mass center to approximate the viscous friction of the fluid through which the motor protein moves, as shown in Figure 2. To truly assess the effects of viscous friction on a rigid body, one should consider drag which depends on its shape and orientation, but here a simple coefficient of viscous friction is used. Further details are given in Supporting Information A.1. 2.3. Conformational Changes and External Forces. Several forces contribute to protein locomotion, including the conformational changes due to ATP hydrolysis. Regardless of their source and application, their resultant can be transformed into equivalent forces at the protein’s joints. At this point in the development of the multibody model, it is not critical to model all of the forces contributing to protein locomotion in detail. The equivalent forces provide a simple means for generating stepping, which allow examination of contact and impact between the heads and substrate. The viscous friction, binding charge, contact, and random forces are excluded from the equivalent forces because they are modeled explicitly. The equivalent forces are modeled as three moments, one acts between the necks, and two act between the heads and necks, see Figure 3. They are loosely associated with the neck linker force and power strokes at each head; however, they may have several sources beyond ATP hydrolysis. These equivalent forces comprise Γconform in eq (1). (23) Baruh, H. Analytical Dynamics, 1st ed.; WCB MacGraw-Hill: 1999. (24) Mather, W. H.; Fox, R. F. Biophys. J. 2006, 91, 2416–2426.
12976 DOI: 10.1021/la901812k
Bowling et al.
Figure 3. Forces emulating conformational changes and external forces.
Figure 4. Binding charges on myosin V.
A combination of constant neck linker24,25 and power stroke forces is used here to produce locomotion. This is done so that the same forces can be applied to the massless and massive models. The motors are activated when both heads have docked, only one power stroke is active at any given time. Both forces are deactivated when the stepping head’s binding site is 5 nm from the actin binding site. Allowing the actin binding site to pull the head in to bind and dock, yields a more robust control, as opposed to steering the head all the way in. In addition, the power stroke is deactivated when the angle between a head and (25) Jamali, Y.; Lohrasebi, A.; Rafii-Tabar, H. Phys. A 2007, 381, 239–254.
Langmuir 2009, 25(22), 12974–12981
Bowling et al.
Article
Figure 6. Brownian motion for myosin V.
Figure 5. Contact forces on myosin V.
its neck is 90° assuming this is the maximum travel of the power stroke. The neck linker force is deactivated if the angle between the necks is greater than 75°. The motors apply a constant torque when activated. Further details are given in Supporting Information A.2. 2.4. Binding Charges. In the literature, there are charge potentials proposed which attempt to model the repulsion required to keep the head from penetrating the actin filament. Since the repulsion forces are addressed using contact forces, herein the binding potential only models the attraction between the binding sites using a simple Coulomb-like charge model. The attractive forces resulting from these charges are shown in Figure 4; a detailed discussion is given in Supporting Information A.3. The interaction of the motor protein with the actin filament involves docking of the heads at the binding sites.26,27 Here, this docking is assumed to completely immobilize the bound head and only occurs when the point binding sites on a head and the actin filament closely align. Modeling the head as an ellipse ensures that the head must achieve a particular configuration before docking is allowable, as would occur for the actual protein. When a head docks, it is assumed that the charges involved sum and are neutralized; the charges at the docking site do not affect the undocked head. This is accomplished by setting the charges involved equal to zero. These charges remain zero until the head detaches, after which only the charge on the binding site is set back to its original value. After it is detached, the head must recharge at some point during its motion in order to achieve the step. The head charge is switched back on when its binding site has passed mid step and is 2/3 of the way to the targeted binding site on the actin filament. The head needs to be closer to the new binding site than to the old since the binding charges are a function of distance. If not, the head will be pulled back to the original binding site. (26) Bier, M. Phys. Rev. Lett. 2003, 91, 148104-1–148104-4. (27) Bier, M. Contemp. Phys. 2005, 46, 41–51.
Langmuir 2009, 25(22), 12974–12981
2.5. Contact Forces. Here it is assumed that contact between the heads and the actin filament occurs in a localized region which can be approximated as point contact. Modeling the heads as ellipses aligns with this assumption. In this work the heads and actin are assumed to be frictionless, for the sake of simplicity, thus the horizontal forces illustrated in Figure 5 do not represent tangential friction. They are used to enforce an immobilization of the head when it docks. See Supporting Information A.4 for further details. 2.6. Brownian Motion. Random forces and moments in the model, representing Brownian motion, are implemented as Gaussian white noise. They act at and about the mass center of each body, as shown in Figure 6. See Supporting Information A.5 for further details.
3. First-Order Model Division of eq (1) by the viscous damping coefficient yields 0≈
:: _ qÞ ΓcccB AðqÞqþbðq, ¼ -Dq_ β β
ΓcccB ¼ Γcontact þ Γcharge þ Γconform þ ΓBrown
ð5Þ
ð6Þ
The argument behind the small mass assumption is that the mass terms divided by β are sufficiently small to be omitted; this is clearly untrue for large accelerations. Regardless, 5 yields a first-order model where velocity is directly proportional to force ΓcccB aA ΓcccB ¼ D -1 q_ ≈ q_ ¼ D -1 β β aA
ð7Þ
assuming D is invertible. This assumption is central to the widely used particle models of Langevin28-30 and Fokker-Planck.31,32
(28) Reif, F. Fundamentals of Statistical and Thermal Physics; McGraw Hill: New York, 1965. (29) Yu, J.; Ha, T.; Schulten, K. Biophys. J. 2006, 91, 2097–2114. (30) Rafii-Tabar, H.; Jamali, Y.; Lohrasebi, A. Phys. A 2007, 381, 239–254. (31) Zeldovich, K. B.; Joanny, J. F.; Prost, J. Eur. Phys. J. E 2005, 17, 155–163. (32) Wang, H.; Elston, T. C. J. Stat. Phys. 2007, 128, 35–76.
DOI: 10.1021/la901812k
12977
Article
Bowling et al.
Note that the scaling of the generalized active forces in (7) implies that the motion is driven by small forces. However, the first-order model has a highest derivative of smaller order than the original model in eq (1), referred to as reduction of order, which yields a singular perturbation problem.17 This is an issue because the first order model can only satisfy an initial position condition; dynamic simulations use initial positions and velocities. In this work, the initial velocity condition is used to encode the effect of impact and contact with the environment. Impact forces are of such short duration that their effect is modeled as an instantaneous change in velocity; herein, contact is modeled as a rapid succession of impacts. The numerical integration is stopped and restarted with updated initial velocities. The impact is lost if the initial velocity condition cannot be satisfied.
4. Second-Order Model Within the framework of perturbation theory, eq (5) is the zeroorder perturbation of eq (1). It represents the protein’s slowest changing dynamics.17 However, it indicates a cancellation of generalized active forces which can be eliminated from the model leaving small forces assumed to drive the motion. This is accomplished using generalized active forces expressed as :: _ qÞ ¼ ΓcccB þ Γfriction ð8Þ AðqÞq þ bðq,
Figure 7. Massless myosin V single step. The trailing head maintains contact with the substrate throughout most of the step as compared to the massive protein, (Head D is green, CPUtime = 113 mins, AbsTol = 10-7, RelTol = 10-8, Δt = 0.01 ms, vav = 1.11 nm/ms).
Separating Γfriction into large and small components yields _ Γfriction ¼ -βða1 Dq_ þ a2 DqÞ ¼ -βa1 Dq_ - β2 Dq_
ð9Þ ð10Þ
where a1 þ a2 = 1 and the value of β2 is determined empirically to correlate simulation and experimental data. Combining eq (7) and eq (10) yields ΓcccB - β2 Dq_ ð11Þ Γfriction ¼ -β1 aA where β1 = a1aA and β2 = βa2. This cancels the large friction forces from eq (11) and from the model. Continuing along this line, ΓcccB in eq (8) is defined as ΓcccB ¼ ðβ1 þ 1Þ
ΓcccB aA
ð12Þ
which is decomposed into large and small components. Substituting eq (11) and eq (12) back into eq (8) yields Γ :: _ qÞ ¼ cccB - β2 Dq_ AðqÞq þ bðq, aA
ð13Þ
where the large forces cancel. It is desirable to retain the relationship between the generalized active forces in eq (7) for small accelerations. Omitting the mass terms in eq (13) yields D -1 ΓcccB ð14Þ q~_ ¼ β2 aA Comparing eq (7) and eq (14), q~_ ¼ q_ if aA = β/β2. Using this relation along with β2 = βa2 and a1 þ a2 = 1 yields β1 ¼
β -β2 β2
Thus, β2 is the only freely chosen parameter in eq (13). 12978 DOI: 10.1021/la901812k
ð15Þ
Figure 8. Massive myosin V single step. The trailing head rolls over and lifts away from the substrate very early in the step as compared to the massless protein, (Head D is green, CPUtime = 15 mins, AbsTol = 10-7, RelTol = 10-8, Δt = 0.01 ms, vav = 1.28 nm/ms).
The parameter β2 is chosen empirically to match experimentally observed motions, but should be small to avoid the original proportionality problem. A small β2 implies that the protein experiences small effective viscous friction forces. In turn, this implies that the system might not be over damped, and may display other second order behaviors such as under or critical damping. Others have suggested a similar conclusion.33-35
5. Examples Hereafter, the first and second order models are referred to as the massless and massive proteins respectively. The dynamics of the massless protein are defined by eq (7). Its velocity is adjusted to keep the protein from penetrating the actin filament if there is contact or impact. The massive protein model is as defined in section 4. In the following sections, deterministic protein locomotion is compared to that obtained in the presence of Brownian motion. Matlab’s ode45.m was used; see Supporting Information B and C. 5.1. Deterministic Stepping. The model forces were adjusted to obtain the experimentally observed velocity of (33) Feynman, R. P. Eng. Sci. 1960, 23, 22–36. (34) Sohlberg, K.; Tuzun, R. E.; Sumpter, B. G.; Noid, D. W. Nanotechnology 1997, 8, 103–111. (35) Kang, J. W.; Hwang, H. J. Nanotechnology 2004, 15, 614–621.
Langmuir 2009, 25(22), 12974–12981
Bowling et al.
Figure 9. Final undocked position of massless protein’s head C from Figure 7. The head binding site must lie within a 0.02 nm W 0.001 nm H box centered atop the actin binding site to be considered as docked.
1 nm/ms. This required β2 = 1 ag/ms, with all other parameters of the two models being equal; increasing viscous damping causes the massive protein to move slower. Recall that for the massless protein β = 108 ag/ms. This difference supports the notion that the effective friction felt by the motor protein is quite small and its locomotion may not be over damped. A single step for the massless and massive proteins is shown in Figure 7 and Figure 8. Similar time snapshots are given in each figure. At the beginning of these simulations both heads are docked. The binding charge at the trailing head is deactivated and the conformational forces are activated to initiate a stepping motion. When the trailing head passes mid step and reaches twothirds of the way to the target binding site, its charge is reactivated and it is drawn toward the new binding site. A key difference between Figure 7 and Figure 8 is that the trailing head of the massless protein remains in contact with the actin filament throughout most of the step. After passing head D at mid step, head C begins to roll over and finally leaves the filament. This is not pure rolling, but a rolling-slipping motion. The trailing head of the massive protein in Figure 8 rolls over much earlier and then drifts upward losing contact with the filament. This behavior is more realistic because of viscous friction and the momentum of the head as it rolls over. This motion results in a type of hand-overhand motion, similar to bipedal walking, as discussed in the literature.36,12 The massive protein will repeat this behavior for several steps. The massless protein has difficulty docking at the binding site. Docking occurs when the head’s binding site falls within a 0.02 nm W 0.001 nm H box which sits atop and is centered on the actin binding site. This small docking zone is used so that the head must achieve a particular configuration in order to dock. This mimics the actual docking of the protein which must align itself with the docking site to make a physical connection. The closeup of the massless protein’s head C in Figure 9 shows the head’s contact point and both binding sites. The massless protein’s head undershoots this region but has difficulty correcting the error. The numerical integration significantly slowed down at this point and was stopped. The head (36) Warshaw, D. M.; Kennedy, G. G.; Work, S. S.; Krementsova, E. B.; Beck, S. Biophys. J. 2005, 88, L30–L32.
Langmuir 2009, 25(22), 12974–12981
Article
reached the binding site at t = 32.34 ms and had not docked at t = 37.13 ms after 4.89 ms had elapsed. The massive protein has no difficulty docking at t = 28.08 ms. The head overshoots the binding site but oscillates until docking is achieved. The oscillations are permissible because of the small amount of effective damping, β2 = 1 ag/ms in eq (13), actually experienced by the motor protein. The massless and massive proteins have average speeds of approximately 1.11 nm/ms and 1.28 nm/ms, which are near myosin V’s experimentally observed average speed of 1 nm/ms. However, the massless protein never actually docked. Recall that achieving a comparable speed for both protein’s required β2 = 1ag/ms; increasing β2 slows the massive protein down. These stepping motions are sensitive to changes in the model parameters such as the specific initial configuration, the magnitude of the binding charges, etc. However, Figures 7 and 8 were obtained using identical model parameters, just different modeling approaches, so they do provide an example of the differences that can occur. 5.2. Random Stepping. It is reasonable to assume that Brownian motion would fix the docking problems of the massless protein. In this section, the same sets of random forces are applied to the massive and massless models used in section 5.1 to investigate this possibility. Brownian motion is implemented using Gaussian white noise. With this type of modeling it is expected that the system should return to its original configuration at some point. Thus the conformational forces from the deterministic stepping are used to give the protein some direction. An adaptive numerical integrator will have difficulty with discontinuous random forces, so they are kept constant during a single time step. This introduces an artificial frequency to the random forces, which is purely a product of the simulation technique. Longer time steps produce an effect like strong random gusts of wind, while short time steps produce high frequency vibrations. Here a time step of 0.01 ms is used for a large effect. Further details are given in Supporting Information A.5. Examples of the resulting motion for the massive and massless proteins are shown in Figures 10 and 11. It is difficult to display the actual jittery motions, so only the tracks of the binding site on each head and the load’s mass center, points CE, DE, and LOADo, are shown overlaid on the docked configurations. Note that some impact events are missed by the integrator, so the head penetrates the filament. This is an issue with the time step of the numerical integration coupled with the behavior of the event function. The first thing to notice is the diffusion of the heads out in every direction, as indicated by the tracks of the binding sites. This motion has been experimentally observed.12,22 This diffusion increases the time for completing a step beyond that in the deterministic case. Another difference between Figures 10 and 11 is that the massive protein attenuates the effect of Brownian motion as evidenced by the more continuous appearance of the tracks in Figure 10. This is a characteristic common to most dynamic systems in that they tend to low pass filter forces. This attenuation should make it easier for the massive protein to dock in the presence of Brownian motion. This conclusion should apply regardless of the model used to emulate Brownian motion. The tracks in Figure 11 appear almost as a cloud of points and are much more discontinuous. The massless protein actually does not complete a step in Figure 11. It has difficulty docking, just as it did in Figure 7. The simulation slowed down considerably and was stopped at t = 109.48 ms. The head binding site reached the actin binding site at t = 104.59 ms and attempted to dock for DOI: 10.1021/la901812k
12979
Article
Bowling et al.
Figure 10. Massive protein and Brownian motion. The protein has no difficulty docking, but is not helped a great deal by Brownian motion (Head D is green, CPUtime = 4 hours, AbsTol = 10-5, RelTol = 10-6, Δt = 0.01 ms, vav = 0.605 nm/ms). Table 1. Statistical Data for Stepping of Massless and Massive Proteinsa massless simulations
14
massive 17
steps completed repeated dwelled
23 2 27
100 32 67 CPU time
max av ( dev min
Figure 11. Massless protein under the influence of Brownian motion. The simulation was stopped because it slowed down considerably while head C attempted to dock; 4.89 ms elapsed in the docking attempt, which is the same delay that was encountered without Brownian motion, (Head D is green, CPUtime = 2 hours, AbsTol = 10-5, RelTol = 10-6, Δt = 0.01 ms, vav = 0.34 nm/ms).
22.18 hours 10.75 ( 10 hours 49.06 min
2.09 hours 43.59 ( 23.86 min 18.86 min
dwell time (ms) max av ( dev min
13.66 1.15 ( 2.6 2.16 10-8
0.034 0.0068 ( 0.0072 0.00004
speed (nm/ms)
4.89 ms before the simulation was stopped. This is the same amount of time that elapsed attempting to dock in the deterministic case. Thus the random forces used in this particular case do not assist the massless protein in docking. In fact, the heads often come near binding sites and are pushed away by random forces, increasing the stepping time. The random forces do not appear to help the massive protein, but act more as a disturbance force to the deterministic motion. The random forces acting on the massless protein in Figure 11 are plotted later in Figure 13 in the Supporting Information. The effect of random forces can be reduced by decreasing the step size, but at some point, they will only produce small amplitude vibrations of the protein. This decreased effect might improve the behavior of the massless protein in docking, but would also decrease the amount of experimentally observed diffusion out in all directions exhibited by the simulation. These are interesting issues to explore, but a full characterization of Brownian motion with respect to the proposed model is beyond the scope of this paper. Yet, since random forces are involved, it is useful to run several cases to obtain the nominal behavior. In general, Brownian motion did not fix the docking problem for the massless protein, and does not truly aid the motion of either protein. The number of runs out of 20 that yielded successful stepping are shown in the first row of Table 1. In all of the simulations a limit, between 0.5 ms and 3 ms, was enforced on how long the protein could take to dock. If this limit was exceeded, the protein’s configuration was 12980 DOI: 10.1021/la901812k
max 1.77 2.04 av ( dev 0.71 ( 0.53 0.99 ( 0.44 min 0.05 0.54 a Matlab’s ode45.m was used with relerr = 10-5, abserr = 10-4.
altered in order to help it dock. Sometimes this helped the massless protein continue stepping, but sometimes not. Simulations with excessive run times were terminated. Table 1 shows that the number of successful steps taken versus those that produced a dwell, where the head slowly approaches the binding site but does not dock, are high for the massless protein. A repeat occurs when the protein attempts to step, but ends up back in its original configuration. The massless protein never finished all four steps. The run times shown in Table 1 are significantly longer for the massless protein because its motion would slow down considerably during docking. The elapsed time for the dwells is shown in Table 1 which shows that the massive protein rarely dwells for long, as expected. The average speed for the massive protein is fairly close to the experimentally observed one of 1 nm/ms, but the massless protein loses speed because of all the dwell times. These data suggest that Brownian motion does not improve the behavior of either protein, but acts as a disturbance to their motion.
6. Discussion This section offers explanations for the behaviors illustrated and discussed in the previous section. Langmuir 2009, 25(22), 12974–12981
Bowling et al.
6.1. Reduction of Order. The massless protein cannot satisfy the initial velocity which is used to transmit the effect of impact forces on the massive protein. The ability to satisfy an initial velocity gives the massive protein a velocity memory such that its previous velocity affects its current state; the massless protein only has a position memory. The massive protein remembers its rotational velocity from the previous time step, and unless some force intercedes, it will continue to rotate causing the head to roll over. 6.2. Momentum. Velocity memory gives the protein a momentum that changes continuously, except during an impact event. This quantity is conserved unless the protein is acted on by external forces. Since contact and impact forces are involved here, momentum may or may not be conserved. However, the forces acting on the protein should equal the rate of change of momentum; this is a statement of Newton’s second law. The massless protein’s lack of velocity memory may or may not yield a continuously changing momentum, and the changes in momentum may not correspond to the external forces. 6.3. Two Stage Numerical Integration. Dynamic simulation involves integrating the rate of change of momentum in order to obtain the protein’s current velocity and position. Thus a single integration of contact forces changes velocity, and a second integration is required for their effect to reach the position. This two-step integration causes some attenuation of external forces. The effect of forces acting on the massless protein is transmitted to the position after a single integration, and thus is not attenuated much. 6.4. Inertial Coupling. The mass matrix in a multibody model is often nondiagonal which means there is coupling between accelerations of different parts of the protein. This coupling encodes a portion of the interparticle interactions in the mass matrix so that the effect of impact and contact forces
Langmuir 2009, 25(22), 12974–12981
Article
acting on one particle is transmitted to other particles. This effect is lost in the first-order model. 6.5. First- versus Second-Order Behavior. The massless protein’s difficulties in docking display classic over damped behavior in that the head slowly converges toward the binding site, but may never actually reach it; however, this does not occur in every simulation. This is the usual behavior predicted by first order models. In contrast, due to small effective viscous damping, the second order model allows for the possibility of overshoot and oscillation in order to achieve docking.
7. Conclusions This paper presents a new analytic approach toward addressing the simulation of multitime-scale, rigid, multibody, motor protein models. This approach yields a second order model which is compared against a first order model obtained using the wellknown small mass assumption. When the average speed and overall motion are matched, the two models predict different dynamic behavior. Although the first order model captures some aspects of the overall motion, the behavior predicted by the second order model appears more realistic. The development of the second order model implies that the effective friction forces, as well as the other external forces, in the motor protein environment may actually be quite small. This allows the motor protein to display the full complement of second order behaviors including under and critical damping, even though the large viscous damping coefficient suggests an over damped system. The proposed approach provides a model which can be used as the foundation upon which more complex motor protein models can be built. Supporting Information Available: Details of the model development and of the numerical integration algorithm are provided. This information is available free of charge via the Internet at http://pubs.acs.org.
DOI: 10.1021/la901812k
12981
