A multiscale modeling approach for biomolecular systems

Abstract

This paper presents a new multiscale molecular dynamic model for investigating the effects of external interactions, such as contact and impact, during stepping and docking of motor proteins and other biomolecular systems. The model retains the mass properties ensuring that the result satisfies Newton’s second law. This idea is presented using a simple particle model to facilitate discussion of the rigid body model; however, the particle model does provide insights into particle dynamics at the nanoscale. The resulting three-dimensional model predicts a significant decrease in the effect of the random forces associated with Brownian motion. This conclusion runs contrary to the widely accepted notion that the motor protein’s movements are primarily the result of thermal effects. This work focuses on the mechanical aspects of protein locomotion; the effect ATP hydrolysis is estimated as internal forces acting on the mechanical model. In addition, the proposed model can be numerically integrated in a reasonable amount of time. Herein, the differences between the motion predicted by the old and new modeling approaches are compared using a simplified model of myosin V.

Appendices

Appendix A: The model

This section discusses in detail the different types of force acting on the model presented in Sect. 2.2. The physical parameters of the myosin V mechanical model are given in Table 2.

Table 2 Physical parameters of the Myosin V mechanical model

1.1 A.1 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 Fig. 9. In order to truly assess the effects of viscous friction on a rigid body, one should consider drag that depends on its shape and orientation, but here a simple coefficient of viscous friction is used. For example, the friction force and moment on body A are

$$ {\mathbf{f}}_{A_o} = - \beta_2 \boldsymbol {v}_{A_o} \quad \mbox{and}\quad {\mathbf{m}}_{A} = - \beta_2 \bar{L}_A^2 {\boldsymbol{\omega}}_{A} $$

(23)

where $\mathbf{f}_{i_{o}}$ and m _i are the force and moment acting on body i, and $\mathbf{v}_{i_{o}}$ and ω _i are the translational and rotational inertial velocities associated with body i. The term $\bar{L}_{i}$ is a characteristic length for rotational viscous friction; the characteristic length for a sphere shape body equals its radius, which is the half length of the whole body [78]. Thus, the characteristic length for other shape bodies is approximately considered as the half length of the whole body. The values of these coefficients are listed in Table 2. The viscous friction forces for each body comprise Γ _friction.

Fig. 9

Viscous friction forces acting on myosin V

1.2 A.2 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 the forces contributing to protein locomotion in detail. The equivalent forces provide a simple means for generating stepping, which allows 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 acting between the necks, and two acting between the heads and necks; see Fig. 10. 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. (15).

Fig. 10

Forces emulating conformational changes and external forces

A combination of constant neck linker [66, 87] 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, and 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 its neck is 90^∘ assuming that 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. The equivalent moments acting on the proteins to achieve the step are listed in Table 3. These moments are separated into large and small components as

$$ {\boldsymbol{\varGamma}}_{\mathrm{conform}} = \mathsf{G}^T \left [ \begin{array}{c} \varUpsilon_R \\ \varUpsilon_S \\ \varUpsilon_P \end{array} \right ] ( \beta_1 + 1) = \mathsf{G}^T {\boldsymbol{\varUpsilon}} ( \beta_1 + 1 ) $$

(24)

where ϒ _i is the moment associated with point i, and G ^T transforms the moments into generalized active forces.

Table 3 Equivalent joint moments

1.3 A.3 Binding charges

In the literature, there are charge potentials proposed that 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 Fig. 11.

Fig. 11

Binding charges on myosin V

The interaction of the motor protein with the actin filament involves docking of the heads at the binding sites [13, 14]. 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 detaching, 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 functions of distance. If not, the head will be pulled back to the original binding site.

The position vectors from the heads’ mass centers to their binding sites are expressed as

$$\begin{aligned} {\mathbf{P}}_{\mathrm{C}_{\mathrm{o}} \mathrm{C}_{\mathrm{E}}} = & L_6 \mathbf{C}_1 + L_7 \mathbf{C}_2 \end{aligned}$$

(25)

$$\begin{aligned} {\mathbf{P}}_{\mathrm{D}_\mathrm{o} \mathrm{D}_\mathrm{E}} = & L_8 \mathbf{D}_1 + L_9 \mathbf{D}_2 \end{aligned}$$

(26)

The force on a head’s binding site $\rm C_{E}$ due to site $\rm B_{1}$, for example, is

$$ {\mathbf{F}}_{\mathrm{C}_\mathrm{E} \mathrm{B}_{1}} = \mathrm{f}_{\mathrm{C}_\mathrm{E} \mathrm{B}_{1}} \frac{\mathbf{P}_{\mathrm{C}_\mathrm{E} \mathrm{B}_{1}} }{\| {\mathbf{P}}_{\mathrm{C_E B_1}}\|} $$

(27)

There are a number of charge potentials investigated in the literature for defining the attractive force $\mathrm{f}_{\mathrm{C_{E} B_{1}}}$. These include ratchet potentials [66] and Coulomb-like potentials [32]. A Coulomb-like charge potential [32] is used here:

$$ V(r) = \frac{ C_{\mathrm{C}_\mathrm{E}} C_{\mathrm{B_1}} e^{k(a-r)} }{ 4 \pi\epsilon_0 \epsilon_r (1 + k a) (r+k_o)} $$

(28)

where ϵ ₀, ϵ _r , k, and a are the permittivity of free space, relative electric permittivity, inverse of the Debye length, and excluded volume radius; the values for these constants are given in Table 2. The term r is the distance between the binding sites, and k _o is a penetration depth used to keep the potential from becoming infinite. The terms $C_{\mathrm{C}_{\mathrm{E}}}$ and $C_{\mathrm{B}_{1}}$ are the charges on the binding sites of head C and B₁ on the actin fiber; these values are given in Table 2. The absolute value of this potential as a function of radial distance is illustrated in Fig. 12.

Fig. 12

Absolute value of the charge potential and the associated force

The force can be found as

$$ {\mathrm{f}}_{\mathrm{C}_\mathrm{E} \mathrm{B}_{1}} = -\frac{dV(r)}{dr} = \frac{C_{\mathrm{C}_\mathrm{E}} C_{\mathrm{B_1}} e^{k(a-r)} }{4 \pi\epsilon_0 \epsilon_r (1 + k a) } \frac{ k(r+k_o) + 1}{(r+k_o)^2} $$

(29)

where $r=\| {\mathbf{P}}_{\mathrm{C}_{\mathrm{E}} \mathrm{B}_{1}} \|$ is the distance between the head and actin binding sites at points C_E and B₁ in Fig. 2. Similar forces can be defined between the other binding sites. These charge forces comprise Γ _charge in (15) and are separated into small and large components:

$$ {\boldsymbol{\varGamma}}_{\mathrm{charge}} = (\beta_1+1) Q \mathbf{f}_{\mathrm{charge}} $$

(30)

where the matrix Q transforms the charges into generalized active forces. The force $\mathrm{f}_{\mathrm{C}_{\mathrm{E}} \mathrm{B}_{1}}$ generated by the charge potential as a function of radial distance is illustrated in Fig. 12.

These charges are electrostatic and only depend on the relative positions of the binding sites. They are conservative forces whose values are always known. Thus, it is necessary to ensure that their computed values and the decomposed values match up. This is accomplished by defining $\overline{\mathrm{f}}_{ij}$:

$$ {\mathbf{f}}_{\mathrm{charge}} = [ \overline{\mathrm{f}}_{\mathrm{C}_\mathrm{E} \mathrm{B}_{1}}\quad \overline{\mathrm{f}}_{\mathrm{C}_\mathrm{E} \mathrm{B}_{2}} \quad\cdots\quad\overline{\mathrm{f}}_{\mathrm{D}_\mathrm{E} \mathit{B}_{5} } \quad \overline {\mathrm{f}}_{\mathrm{D}_\mathrm{E} \mathrm{B}_{6}} ]^T $$

(31)

where the $\overline{\mathrm{f}}_{ij}$ are solved such that, for example,

$$ {\mathrm{f}}_{\mathrm{C_E B_1}} = (\beta_1 + 1) \overline { \mathrm{f}}_{\mathrm{C_E B_1}} $$

(32)

and $\mathrm{f}_{\mathrm{C_{E} B_{1}}}$ is determined from (29).

1.4 A.4 Contact forces

1.4.1 A.4.1 Contact points

Here it is assumed that contact between the heads and the actin filament occurs in a localized region that can be approximated as point contact. Modeling the heads as ellipsoids 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 Fig. 13 do not represent tangential friction. They are used to enforce an immobilization of the head when it docks.

Fig. 13

Contact forces on myosin V

As stated, the heads are modeled as ellipsoids:

$$ \frac{ x_C^2 }{ r_{\mathrm{head}_1}^2 } + \frac{ y_C^2 }{r_{\mathrm{head}_2}^2 } + \frac{ z_C^2 }{ r_{\mathrm{head}_3}^2 }= 1 $$

(33)

where the position vector between point C_o and any point C_C on the surface of the head is expressed as

$$ {\mathbf{P}}_{C_o C_C} = x_C \mathbf{C}_1 + y_C \mathbf{C}_2 + z_C \mathbf{C}_3 $$

(34)

and likewise for head D. Table 2 gives the values of $r_{\mathrm{head}_{1}}$, $r_{\mathrm{head}_{2}}$, and $r_{\mathrm{head}_{3}}$ assuming that both heads are identical.

At each step in the simulation, the contact points between each head and the actin filament are determined. The approach used to determine the contact forces at these points is presented in [20] and will not be duplicated fully here. Tangential forces associated with sliding or sticking friction are omitted.

1.4.2 A.4.2 Complementarity conditions

The complementarity conditions are a well-known means for describing the relationships between friction, contact forces, velocities, and accelerations [2]. Assuming that the distance between the impacting points equals zero, the conditions are dependent on the value of the preimpact normal velocity and acceleration as

$$ \left \{ \begin{array}{l@{\quad}c@{\quad}l} v_{n_i}(t) < 0 & & \mathrm{impact or contact} \\ v_{n_i}(t) = 0 & \mathrm{and} & \left \{ \begin{array}{l@{\quad}l} \dot{v}_{n_i}(t) \leq0 & \mathrm{contact} \\ \dot{v}_{n_i}(t) > 0 & \mathrm{separation} \end{array} \right . \\ v_{n_i}(t) > 0 & & \mathrm{separation} \end{array} \right . $$

(35)

A transition occurs when the preimpact normal velocity equals zero. The preimpact acceleration must be checked to determine whether impact forces will exist.

When the pre-impact normal velocity indicates impact or contact, its postimpact value is

$$ v_{n_i}(t+\epsilon) = -e_{n_i} v_{n_i}(t) $$

(36)

Herein the nonpenetration condition is modeled by $e_{n_{i}} = 0$.

1.4.3 A.4.3 Contact force calculation

In brief, the method treats contact as a succession of small impacts. Impact is treated as a discrete event with an instantaneous velocity change. At each step in the numerical integration, an impact law is used to determine an instantaneous change in velocity. This law uses a coefficient of restitution, which is set to zero for the motor protein case. This means that the head does not rebound but can remain attached to the substrate after impact. This approach also involves the calculation of the contact and impact forces that comprise Γ _contact. These forces are necessary, in addition to the velocity resolution, in order to keep the head from penetrating the actin filament. See [20] for further details.

The contact forces are shown in Fig. 13 and defined as

$$\begin{aligned} \boldsymbol {\varGamma }_{\mathrm{contact}} = J_p^T \mathbf{F}_{\mathrm{contact}} = J_p^T \left [ \begin{array}{c} \mathbf{f}_{C_{C}} \\ \mathbf{m}_{C} \\ \mathbf{f}_{D_{C}} \\ \mathbf{m}_{D} \\ \mathbf{f}_{R} \end{array} \right ] \end{aligned}$$

(37)

$$\begin{aligned} \boldsymbol {v}_p = \left [ \begin{array}{c} \boldsymbol {v}_{C_{C}} \\ \boldsymbol {v}_{D_{C}} \\ \boldsymbol {v}_{R} \end{array} \right ] = J_p \ \dot{\mathbf{q}} \end{aligned}$$

(38)

where v _p contains the inertial velocities of the contact points C _C and D _C , as well as that of point R. Although all possible contacts are shown in (38), the ones that actually contact/impact are selected using the impact Jacobian J _p .

The contact forces are calculated as

$$\begin{aligned} {\mathbf{F}}_{\mathrm{contact}} =& - \bigl( J_p A^{-1} J_p^T \bigr)^{-1} \bigl( \dot{J}_p \dot{\mathbf{q}} + J_p A^{-1} ( \boldsymbol {\varGamma }_{\mathit{fccB}} - \mathbf{b} ) \bigr) \end{aligned}$$

(39)

$$\begin{aligned} \boldsymbol {\varGamma }_{\mathit{fccB}} =& \boldsymbol {\varGamma }_{\mathrm{friction}} + \boldsymbol {\varGamma }_{\mathrm{charge}} + \boldsymbol {\varGamma }_{\mathrm{conform}} + \boldsymbol {\varGamma }_{\mathrm{Brown}} \end{aligned}$$

(40)

Since F _contact depends on the other forces, its large portion is found by substituting the other large forces into (39) and solving.

1.5 A.5 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 Fig. 14.

Fig. 14

Brownian motion for myosin V

The random forces and moments shown in Fig. 14 representing Brownian motion are defined, for example, as

$$ \begin{aligned} {\mathbf{f}}_{Co} & = C_{o1}(t) \mathbf{N}_1 + C_{o2}(t) \mathbf{N}_2 + C_{o3}(t) \mathbf{N}_3 \\ {\mathbf{m}}_C & = \bar{L}_C C_{o1}(t) \mathbf{N}_1 + \bar{L}_C C_{o2}(t) \mathbf{N}_2 + \bar{L}_C C_{o3}(t) \mathbf{N}_3 \end{aligned} $$

(41)

where $\bar{L}_{C}$ is a characteristic length of body C. Similar forces act on the other bodies. The C _oi (t) represent forces produced by randomly fluctuating thermal noise. Each component of the random force and moment is treated independently as a normally distributed random variable [66]. They have the following expectations E[⋅] or weighted average values:

$$ E \bigl[ C_{oi}(t) \bigr] = \bigl\langle C_{oi}(t) \bigr\rangle = 0 = \mu $$

(42)

and are governed by a fluctuation–dissipation relation expressed as

$$ E \bigl[ C_{oi}(t_1) C_{oj}(t_2) \bigr] = 2 \beta k_B \mathrm{T} \delta(t_1-t_2) \delta_{i,j} $$

(43)

where k _B and T are the Boltzmann constant and absolute temperature [59, 66]. The relation in (43) implies that there is no time dependency between the random process over time; the random sequence of forces does not repeat regularly.

In addition, (42) and (43) imply

$$ E\bigl[ C_{oi}^2(t) \bigr] = 2 \beta k_B \mathrm{T} = \mathrm {Var}\bigl(C_{oi}(t)\bigr) = \sigma^2 $$

(44)

which is the variance of C _oi . Thus, the C _oi can be generated using the Matlab function normrnd(μ,σ,…), which generates random variables with normal distribution.

The collection of random forces comprise Γ _Brown. These randomly fluctuating discontinuous functions slow numerical integration so that each random variable is held constant during a single integration step; the random variable is updated at the beginning of each step. Thus, the value of each random variable is known before the integration step, and the decomposed value of the random force must equal it. This is accomplished by defining

$$ \boldsymbol {\varGamma }_{\mathrm{Brown}} = (\beta_1 +1) {R}_{nd} \mathbf{r}_{nd} = (\beta_1 +1) {R}_{nd} \left [ \begin{array}{c} \bar{C}_{o1} \\ \bar{C}_{o2} \\ \bar {C}_{o3} \\ \vdots \end{array} \right ] $$

(45)

where R _nd transforms the random forces into generalized active forces, and

$$ C_{oi} = (\beta_1+1) \bar{C}_{oi} $$

(46)

and likewise for the other random forces. An example of the random forces used is given in Fig. 15.

Fig. 15

Random forces acting on head C in Fig. 8

Appendix B: Simulation algorithm

Here it is assumed that the equations of motion in (15) have been expressed as a system of first-order differential equations in order to integrate them. This can be accomplished using a simple transformation such as

$$ \dot{\mathbf{q}} = \mathbf{u} \qquad \ddot{\mathbf{q}} = \dot{ \mathbf{u}} $$

(47)

such that

$$ \dot{\mathbf{x}}(t) = \left [ \begin{array}{c} \dot{\mathbf{q}}(t) \\ \dot{\mathbf{u}}(t) \end{array} \right ] = \left [ \begin{array}{c} \mathbf{u}(t) \\ A^{-1} ( {\sum}{\boldsymbol{\varGamma}} - \mathbf{b} ) \end{array} \right ] . $$

(48)

Nominally at each time step, $\dot{\mathbf{x}}(t)$ is calculated and used to increment x(t) yielding x(t+Δt). This calculation is modified here to enforce nonpenetration and docking constraints. At each time step, if a contact or impact is identified, the postimpact state $\dot{\mathbf{x}}(t+\epsilon)$ is calculated, and the numerical integration is carried on from there; ϵ is considered an infinitesimally small amount of time.

The impact/contact resolution process involves identifying when an impact or contact occurs, resolving the postimpact/contact velocity, and possibly restarting the integration with the updated postimpact velocity. The details involved in addressing nonpenetration include: (1) treating contact as a succession of impacts, (2) using an impact law that addresses simultaneous impacts and contacts, (3) using an adaptive and event-driven numerical integration scheme, and (4) calculating contact forces. The details of this process, discussed in another paper [20], are too lengthy to duplicate here, so only a brief description of them is given below. The parameters for the simulation are given in Table 2.

2.1 B.1 Contact as a succession of impacts

Herein, a distinction is drawn between impact and contact forces, which are of short and long duration, relative to each other and to the time step of the numerical integration. Impact forces correspond to abrupt changes in velocity, as occurs during docking with the actin filament, whereas contact forces do not. Treating contact as a succession of impacts allows them both to be addressed within the same mathematical framework. This facilitates the development of an impact law, discussed in detail in another paper [20], which can address simultaneous contact and impact, as occurs during protein locomotion.

2.2 B.2 Simultaneous impact law

This law is intended to enforce nonpenetration and other motion constraints. It is assumed, for simplicity, that the head and actin filament are lubricated sufficiently by the intracellular fluid to make them frictionless. When undocked, only nonpenetration constraints are enforced on the heads. When docked, a head is immobilized by enforcing nonpenetration, nonsliding, and nonrotation constraints. Point R is prevented from penetrating the filament as well. Again, for simplicity, it is assumed that when a head impacts the actin filament, it does not rebound. The associated constraint forces are calculated for use in determining contact and impact events.

2.3 B.3 Adaptive event-driven simulation

An adaptive numerical integrator is used because it can address phenomena occurring at large and small time scales while maintaining a specified accuracy. Matlab’s ode45.m is used here along with an event function for finding contact and impact events. The simulation is stopped and restarted at each impact event. However, since contact is treated as a succession of impacts, this would yield a very slow numerical integration.

This is avoided by monitoring contact forces in the event function. It is clear when an impact event occurs because the contact forces jump from zero to nonzero values; during contact, they sustain nonzero values. The discontinuities in force and velocity caused by an impact event will create difficulties for an adaptive numerical integrator, so it must be stopped and restarted. During contact, the forces and velocities are more continuous and do not create large errors in the numerical integration, and thus stopping and restarting the simulation is not required. This practice increases the speed of the simulation considerably.

Appendix C: Error analysis

This section discusses in detail the concept of error analysis used in the numerical integration. It is difficult to assess the error in the first-order model because it inherently involves abrupt, discontinuous changes in forces and velocities. This creates difficulties in numerically integrating the checking functions used herein, which increases the run time significantly. Therefore, this error analysis is only applied to the second-order model. The correctness of the second-order model should carry over to the first-order model since they are identical, except that one omits the generalized inertia forces.

Since the motor protein’s environment contains nonconservative forces, such as dissipative viscous friction, the system energy is not conserved. Thus the work–energy theorem [10] is used to account for the effect of nonconservative forces:

$$ T_1 + W_{1 \rightarrow2} = T_2 $$

(49)

where T _i is the kinetic energy at state i, W _1→2 is the work done in moving from state 1 to 2 including the work done by nonconservative forces. The total work can be calculated as

$$ W_{1 \rightarrow2} = {\displaystyle\int^2_1} \dot{W}_{1 \rightarrow 2} \, dt , \qquad \dot{W}_{1 \rightarrow2} = {\textstyle\sum} {\mathbf{F}} \cdot \dot{\mathbf{q}} $$

(50)

where ∑F represents the generalized active forces in (15). The time derivative $\dot{W}_{1 \rightarrow2}$ can be numerically integrated along with the model. Subtracting this work from the kinetic energy at state 2 should yield a constant number to within the specified error tolerances. In order to examine the errors, a checking function is defined as follows:

$$ {\mathrm{check}}_{1\rightarrow i} = T_i - ( T_1 + W_{1 \rightarrow i} ) $$

(51)

The evaluation of the checking functions involves a difficulty due to the calculation of the work done by the contact forces during the stepping. This is because of the discontinuities that contact/impact create in velocities and forces. The work done by contact forces cannot be captured in (50) while the change in velocities appears in the kinetic energy; that leads to a gap (jump) in the checking function (51). In order to address this issue, the difference between the kinetic energy at the moment of contact/impact, which equals to the nonconservative work, is added to the check function when the numerical is stopped and restarted at the impact event, see Appendix B, to compensate for the gap.

The check function is examined for a single deterministic step, corresponding to Fig. 5, of motor protein locomotion as shown in Fig. 16. The deterministic step is chosen because the simulation of Brownian motion is very discontinuous and difficult to integrate along with the check function. However, the deterministic step gives confidence in the results from the simulations with Brownian motion. It involves large numbers, which may give the appearance that the integrator is not maintaining the desired accuracy, even though it may be for the variables in $\dot{\mathbf {x}}$; see (48). This situation can be remedied by normalizing the checking functions using a characteristic system value. The system energy E _i at state i is defined as

$$ E_i = T_i + V_i $$

(52)

where V _i is the potential energy. Since energies are being evaluated, as shown in Fig. 17, the system energy E _i is used for this purpose. The characteristic energy is taken from Fig. 17 as 10⁵.

$$ {\mathrm{chk}}_{1 \rightarrow i} = \frac{\mathrm{check}_{1 \rightarrow i}}{E_i} $$

(53)

Fig. 16

The checking function of myosin V’s locomotion during one deterministic step that includes drifting toward actin filament, docking to binding sites, and taking one step

Fig. 17

Total energy of myosin V during one step

The normalized checking function chk_1→i is shown in Fig. 18. Notice that the scaled checking function satisfies the specified AbsTol and RelTol, whereas the unscaled one, Fig. 16, does not. The check function in Fig. 18 shows that the desired tolerance is maintained even in the presence of contact and impact forces. This gives confidence in the model’s correctness and the techniques involved in the numerical integration.

Fig. 18

The normalized checking function remains beneath the specified AbsTol=10⁻⁵ and RelTol=10⁻⁶ during the simulation

Appendix D: Units

This appendix gives a brief summary of the units used in the simulation. The unit basis consists of the attogram (ag), the nanometer (nm), the millisecond (ms), and the radian (rad). With this choice the force unit becomes

$$ 1 \mathrm{N} = 1~\frac{\mathrm{kg}\cdot \mathrm{m}}{\mathrm{s}^2} \ \frac{10^{21}\,\mathrm{ag} }{\mathrm{kg}} \ \frac{ 10^9 \,\mathrm{nm} }{\mathrm{m }} \ \frac{\mathrm{s}^2 }{ (10^3~\mathrm{ms} )^2 } = 10^{24}~\frac{\mathrm{ag} \cdot \mathrm{nm} }{\mathrm{ms}^2} $$

(54)

and the unit of work becomes

$$ 1~\mathrm{J} = 1\ \frac{\mathrm{kg \cdot m}^2}{\mathrm{s}^2} = 10^{33}~\frac{\mathrm{ag} \cdot \mathrm{nm}^2 }{\mathrm{ms}^2} $$

(55)

These values are chosen so that the masses have values near unity and to improve the stability of the numerical integration.