Computer simulation of flagellar movement X: Doublet pair splitting and bend propagation modeled using stochastic dynein kinetics

Experimental observations on cyclic splitting and bending by a flagellar doublet pair are modeled using forces obtained from a model for dynein mechanochemistry, based on ideas introduced by Andrew Huxley and Terrill Hill and extended previously for modeling flagellar movements. The new feature is elastic attachment of dynein to the A doublet, which allows movement perpendicular to the A doublet and provides adhesive force that can strain attached dyneins. This additional strain influences the kinetics of dynein attachment and detachment. Computations using this dynein model demonstrate that very simple and realistic ideas about dynein mechanochemistry are sufficient for explaining the separation and reattachment seen experimentally with flagellar doublet pairs. Additional simulations were performed after adding a “super‐adhesion” elasticity. This elastic component is intended to mimic interdoublet connections, normally present in an intact axoneme, that would prevent visible splitting but allow sufficient separation to cause dynein detachment and cessation of shear force generation. This is the situation envisioned by Lindemann's “geometric clutch” hypothesis for control of dynein function in flagella and cilia. The simulations show abrupt disengagement of the “clutch” at one end of a bend, and abrupt reengagement of the “clutch” at the other end of a bend, ensuring that active sliding is only operating where it will cause bend propagation from base to tip. © 2014 The Authors. Cytoskeleton, Published by Wiley Periodicals, Inc.


Introduction
M ovements of cilia and flagella are generated by dynein motor enzymes, operating within a cytoskeleton known as the axoneme. The axoneme is a cylindrical array of parallel elements referred to as outer doublets. This work uses computer simulations to examine the function of structural components that hold the axonemal doublets together. In a working cilium or flagellum, these components maintain a center-to-center spacing between adjacent doublets that is close to the resting value of about 60 nm, while allowing sliding between adjacent doublets that can be more than 6 100 nm. Early work described experimental conditions where a portion of an axoneme can transiently separate into two or three bundles of doublets, which then reassociate and continue to generate active sliding [Sale, 1986;Brokaw, 1986a]. Kamiya and Okagaki [1986] and Aoyama and Kamiya [2005] used enzymatic digestion to partially deconstruct Chlamydomonas axonemes, and observed pairs of doublets that remained firmly connected at their basal ends. Shear forces and sliding generated by dyneins on one doublet caused cycles of bending, separation, and reassociation. These observations demonstrated a role for dyneins in holding doublets together and also demonstrated that inactivation of dynein shear force by separation of the doublets operates as a built-in "control mechanism" that can produce oscillating movements.
Computer simulations [Brokaw, 2009] verified these interpretations by reproducing the cyclic splitting and reassociation of a doublet pair, using simple mathematical models for dynein shear forces and adhesive forces. The simulation programming is now extended to a more realistic model, which calculates dynein attachment and detachment from chemical kinetics and allows bending of both doublets. The ability of this computer modeling to simulate the experimental results demonstrates that both structural and motor functions of dynein can be explained by conventional ideas about the chemical kinetics of dynein motors. In the two-doublet experimental system [Aoyama and Kamiya, 2005], some structural components were eliminated by the enzymatic digestion used to deconstruct the axoneme. The current simulation programming has also been expanded to mathematically replace some functions of the missing components, so that it is possible to examine the operation of a "geometric clutch" [Lindemann, 1994a] during bend propagation.
In the doublet pairs observed experimentally, dyneins on the A tubule of one doublet interact transiently with substrate sites on the B tubule of the other doublet, so these two doublets will be identified here as the A and B doublets. In the computer model, dynein motors are located at 24 nm intervals in two rows along the A doublet, with a 4-nm longitudinal offset between the rows. The mathematical assumptions used for the motor are pictured in Fig. 1. A motor can attach transiently to one of the substrate sites on the B doublet, which are spaced at 8 nm intervals along the length. This attachment is mediated by the motor stalk, which transmits motor force to the B doublet. The model for dynein motor function is similar to a model used successfully for computer simulations of flagellar movement [Brokaw, 1999], except that the stalk has a fixed length (10 mm) and is allowed to pivot at its attachment points on the motor and on the B tubule, through an angle of up to 6 0.45 rad from a perpendicular to the A doublet. The "power stroke" of the model is a 9 nm movement of the A end of the stalk within the motor, parallel to the A doublet, toward the distal end of the doublet. When this movement is restrained, the motor exerts a shear force determined by the distance, x, from the end stroke position and the elastic force constant, FKX, of the motor. An additional feature introduced here is compliance of the motor in the y direction, perpendicular to the A doublet. This compliance has an elastic constant, FKY. The y movement is unstrained when the attachment point for the A end of the stalk is 15 nm from the surface of the A doublet. If this distance is greater than 15 nm, there is an adhesive force that can act through the stalk to decrease the separation between the doublets. If this distance is less than 15 nm, there is a compression force that can push the doublets apart. A nonlinear expression, detailed in Methods, is used for adhesive force at distances greater than 16.5 nm.
Modeling of dynein force production uses a 5 state ATPase cycle (Fig. 1). Shear force and movement are produced by a power stroke associated with the transition from weakly binding state 4 to strongly binding state 5. Further details are presented in Methods, where parameters required to complete the dynein model are given in Table I.

Steady-State Behavior of the Dynein Model
Computer simulations of the dynein model at 0.5 mM ATP were performed with various viscous loads, to obtain force and velocity under steady-state conditions (Fig. 2). Typical force-velocity behavior, as seen in computations with the original 2-state Huxley [1957] model for myosinactin interaction in skeletal muscle, is shown in Fig. 2A. The parameters were chosen so that the sliding velocity at 0 load was 18.5 mmÁs 21 at 0.5 mM ATP, and the K m for decrease with ATP concentration, 0.18 mM, matched the values measured [Kurimoto and Kamiya, 1991] for unloaded sliding disintegration under the same conditions used for the doublet pair experiments. Parameter choices that give an average force of 2.7 pN/motor at 0 velocity give realistic behavior with doublet pair simulations. Force The 5-state ATPase cycle for the dynein model. This cycle is also shown in Fig. 3 of Roberts et al. [2013], with structural diagrams for a cytoplasmic dynein motor. Unshaded states 2 and 3 are detached states, state 4 is the weakly attached state, and states 5 and 1 are the strongly attached, force-producing states. The cycle follows the heavy arrows, except when mechanical detachment occurs from the strongly bound states, through states 6 and/or 7. ATP binding is required for the 1->2 and 7->2 transitions. ATP hydrolysis accompanies the 2->3 transition. ADP release accompanies the 5->1 and 6->7 transitions. The diagrams in the lower part of the figure provide a visual interpretation of the mathematical description of motor function, for the detached states 2 and 3. The green sector represents the range of stalk tilt. The motor itself is colored red, with the A end of the stalk shown in the two positions corresponding to prepowerstroke and postpowerstroke conformations when the motor is unstrained (detached). The blue stack represents the compliant attachment of the motor to the A doublet, allowing up and down motion (the y direction) without lateral motion or tipping. When a motor is weakly attached to a site on the B doublet in state 4 and transits to strongly attached state 5, it also transits to the postpowerstroke conformation and drags the B doublet to the right, toward the distal end of the A doublet. at 0 velocity increased gradually to 3.0 pN/motor at 0.05 mM ATP.
During unloaded sliding computations, with no forces applied in the y direction, the model maintains a doublet surface separation close to 24 nm. This is the expected separation if the stalks of attached dyneins are at close to their maximum tilt of 0.45 rad. When the surface separation is held at a greater distance, the dyneins are strained in the y direction, which increases detachment rates from state 4 to state 3, state 5 to state 6, and state 1 to state 7. Figure 2B shows results from computations of forces at 0 sliding velocity for a range of values of doublet surface separation. A separation of 29 nm, 5 nm beyond the normal operating distance, is sufficient to reduce dynein attachment to a level that effectively eliminates dynein force production.

Cyclic Sliding and Splitting of a Doublet Pair
A typical result from doublet pair computations using this dynein model is shown in Fig. 3A and Supporting Informa-tion Movie M1. The series of images covers one cycle, beginning just before separation of the distal end of the doublet pair (image A1). After complete separation, the doublet pair enters an association phase (A phase). In the A phase, the doublets relax toward elastic equilibrium, which brings the doublets closer together near the basal end of their separation. This proximity favors reattachment of dyneins to sites on the B doublet, as propagation of an association transition (images A2 and A3). These attached dyneins then develop the shear force that causes buckling of the A doublet, to produce the small separation visible in image A4. This event begins the short propagation phase (P phase), in which a propagating association transition is followed by a propagating dissociation transition. The P phase continues until the association transition reaches the distal end, beginning the dissociation phase (D phase, images A5 to A11). These basic features are the same as in Fig. 1 of the previous modeling paper [Brokaw, 2009] and Fig. 2 of the experimental paper [Aoyama and Kamiya, 2005]. Dyneins remaining attached at a separation of 35 nm were less than 1 per cycle; these were programmatically detached.
As a separated region propagates, during the P and D phases, it elongates and increases in amplitude. These changes are associated with active sliding in the distal associated region. The total shear measured at the tip, at end of  the D phase, was 1.55 mm in the cycle shown in Fig. 3A, and averages 1.48 mm (42 cycles, SD 5 0.08 mm). Part of the shear at the tip results from attachment of doublet A to a bent doublet B, which reaches an angle of about 1.5 rad at the distal end. Subtracting this shear (0.09 mm) from the average total shear at the distal end leaves a shear of 1.4 mm, which represents the sliding that produces the extra length of the A doublet curve in the separated region. All results given below for total shear are the raw values, without adjustment for doublet B curvature. A consistent result with these models is that the active sliding in the associated region exhibits a very constant velocity (shown in Fig. 3 of Brokaw [2009]). In the example in Fig. 3A, the sliding velocity was 17.9 mmÁs 21 until the total sliding reached 1.45 mm, and then decreased noticeably at the end of the D phase. The closeness of this sliding velocity value to the unloaded sliding velocity of 18.5 mmÁs 21 ( Fig. 2A) indicates that only a very low shear force needs to be generated by active sliding during the D phase. This is consistent with computations (not shown) that show no significant increase in the elastic energy stored in the bent doublets as the separation enlarges and propagates. Little or no shear force is required to increase the size of the bend of doublet A, because its curvature decreases as it elongates.
Propagation of association and dissociation transitions is shown more quantitatively in Fig. 4, which can be compared with Fig. 4 of the experimental paper [Aoyama and Kamiya, 2005]. The transition criterion used is a doublet separation of 40 nm, which is 11 nm beyond the point at which dynein shear force falls to 0 (Fig. 2B). Beyond 30 nm, the separation rises to 100 nm in less than 1 mm (not shown). The separation that is required for visualization of transition positions in the experimental photographs is not known. Longer records (not shown) reveal considerable variation from one cycle to another, consistent with the stochastic behavior of the dynein motors. The association transition begins to propagate before the end of the D phase. It accelerates after the end of the D phase, to an apparently constant velocity as the transition propagates from 2.5 to 9 mm, followed by a further acceleration that takes it to the distal end. As shown in Fig. 4, association transition velocities were obtained from a linear fit to the time series of transition positions between 2.5 and 9 mm from the basal end. For the cycle in Fig. 3A, a velocity of 396 mmÁs 21 was measured. The average for 35 cycles at 0.5 mM ATP with this model was 404 mmÁs 21 (SD 5 24), essentially the same as the experimental value [Aoyama and Kamiya, 2005] at 0.5 mM ATP.
The dissociation transition propagates at a lower and gradually decreasing velocity. In some cycles, such as the one shown in Figs. 3A and 4, there appears to be a clear transition between a relatively constant fast velocity during the propagation phase, and a slower propagation velocity after the end of the propagation phase. In many other cycles, such as the second cycle in Fig. 4, the transition appears to be more gradual, and obtaining a velocity by a linear fit to points in the P phase is problematic. This group of points was extended, or occasionally reduced, if that improved the quality of the linear fits, as indicated by the R value of the fitting process, but some curves did not provide an acceptable fit. The linear fit in Fig. 4 gave a velocity of 199 mmÁs 21 for the early fast velocity for the cycle in Fig.  3A; the average for this model at 0.5 mM ATP was 214 mmÁs 21 (23 cycles; SD 5 27). This value is intended to be comparable with the value of about 200 mmÁs 21 given by Aoyama and Kamiya [2005] for the early, fast phase of dissociation propagation.
Between about 4.5 and 8 mm from the basal end the velocity of the dissociation transition usually appears sufficiently constant to allow a linear fit to provide a characteristic velocity. A velocity of 78 mmÁs 21 is obtained for the example shown in Figs. 3A and 4. The mean velocity for this model was 85 mmÁs 21 (31 cycles; SD 5 7). The propagation velocity becomes erratic in the final 2 mm of length near the distal end. The experimental results show that propagation of the dissociation transition speeds up as it reaches the distal end. To reproduce this, the model is The vertical lines at the top delineate the association, propagation, and dissociation phases, indicated by A, P, and D. Association transition positions are shown by small black dots, and larger black dots for positions between 2.5 and 9 mm, which were used for linear fitting shown by the black line. A propagation velocity of 396 mmÁs 21 was obtained for the association transition. The dissociation transition positions are also shown by small dots, with larger red or green dots indicating points used for linear fits. The red points, during the P phase, were used to obtain a propagation velocity (199 mmÁs 21 ) for the early, fast period of dissociation propagation. The green dots, in the D phase between positions 4.5 and 8 mm, were used to obtain a propagation velocity (78 mmÁs 21 ) for the intermediate period of dissociation transition propagation. Time was measured from the beginning of a data recording, not all of which is shown here. modified by eliminating dynein motors from one of the rows on the A doublet in the last 20% of the length. This makes the behavior near the distal end more consistent, usually with an increase in velocity, as shown in the first cycle in Fig. 4.
For comparison, results are shown in Fig. 3B and Supporting Information Movie M2 for a model that is identical, except for omission of the bending moment that acts on the doublet pair as a result of shear forces between the doublets. In this reduced model, as in the Brokaw [2009] model, the only moments are those created within each doublet, when longitudinal force acts along a bent doublet. The obvious difference between these two versions is the reduced bend of doublet B in Fig. 3B. The presence of doublet pair bending moment is also revealed in Fig. 3A when the doublets are associated near the basal end, with minimal sliding and near maximal shear force. In image A3 in Fig.  3A, a concentrated bend is visible near the base. This bend increases the moment (or t-force) that initiates buckling and separation of the doublets. Separation is distinctly visible in image A4, but barely visible in image B4, even though the length of the associated, force-producing, region is similar in both cases. In other respects, the patterns are remarkably similar. The propagation velocity in the early, fast phase of dissociation transition propagation decreased less rapidly than the unloaded sliding velocity and a double reciprocal plot indicates a K m (ATP concentration for half-maximal velocity) of 0.10 mM, slightly higher than the experimental value of 0.07 mM measured by Aoyama and Kamiya [2005].

Results at Reduced ATP Concentration
The propagation velocity in the intermediate phase of dissociation transition propagation decreased at lower ATP concentration more rapidly than the unloaded sliding velocity. A double reciprocal plot indicates a K m of 0.37 mM, approximately twice the value of 0.18 mM for the unloaded sliding velocity of the dynein motors. As a consequence, at lower ATP concentrations there is more time in the cycle for sliding to produce a larger final shear. The final shear increases from 1.5 to 2.15 mm as the ATP concentration is reduced. This increase is likely to cause a decrease in association transition velocity, because the final shear leaves the doublets farther apart at the beginning of the association phase. No experimental data are available for comparison with final shear results from the simulations.

Results With Different Values for the Adhesive Force Constant
A value of 0.75 pNÁnm 21 was used for the adhesive force constant, FKY, for the results shown in Figs. 3-5. Figure 6 shows results from doublet pair computations at 0.5 mM ATP with values for the adhesive force constant, FKY, from 0.2 to 2.0 pNÁnm 21 . The association transition velocity increases substantially as FKY is increased, suggesting that increased adhesive force facilitates association, by newly attached dyneins pulling the B doublet closer to the A doublet. There is also a modest decrease in dissociation transition velocity as FKY is increased. Both these effects suggest that the role of FKY in keeping the doublet separation close to normal is more important than the increased staininduced detachment rate at larger separations. Figure 6 shows this decrease in velocity for the propagation velocity in the intermediate phase of dissociation. This decrease lengthens the D phase and decreases the cycle frequency (not shown). With more time for sliding to occur, and a constant sliding velocity, the final shear displacement increases from 1.32 mm at FKY 5 0.2 to 1.58 mm at FKY 5 2.0. The product of dissociation velocity and final shear displacement, shown by the points on the dashed line, is essentially constant in this range, demonstrating an inverse relationship between dissociation propagation velocity and final shear displacement.
Further computations used an additional modification of the nonlinear formulation used for the adhesive force, so that the separation at which the forces become 0 remained close to 29 nm, as in Fig. 2B, over the full range of FKY values. The results were not significantly different from those in Fig. 6, except for a slightly higher value (318 mmÁs 21 ) for the association transition velocity at FKY 5 0.2 pNÁnm 21 .

Additional Results for Cyclic Sliding and Splitting Not Presented in Detail
A modified dynein model having 0 stalk tilt angle required small increases in the rate constants for state 5 to state 1 to state 2 to maintain the maximum steady state sliding velocity at 18.5 mmÁs 21 , as for the model in Fig. 2. The force at 0 velocity increased to 3.2 pN/motor. Additional modifications were required to obtain doublet pair computation results similar to Fig. 3: the adhesive force constant, FKY, was increased to 2.0 pNÁnm 21 and the elastic bending resistance of each doublet was decreased by 20%. This model had lower early dissociation velocities (214%) but the other quantitive parameters differed from the results for the model in Fig. 3A by 5% or less. Similar results were obtained with an intermediate stalk tilt model, possibly more realistic, with motor stalk tilt limited to 6 0.20 rad.
Doublet pair computations were performed with dynein models with modified values for detachment from the strongly attached states 5 and 1 to states 6 or 7. Reducing the 0-strain rate constants for detachment from these states required only small adjustments in the rate constants for state 5 to state 1 to state 2 to maintain steady state performance similar to Fig. 2A. Even with these detachment rate constants reduced from 2 s 21 to as low as 0.01 s 21 , slightly lower than measured values for myosin-actin detachment [Nishizaka et al., 2000], there were only small quantitative changes: 10% or less for total shear, association velocity, and intermediate dissociation velocity, and 221% for initial dissociation velocity. Changes of similar magnitude were obtained with these rates increased from 2 s 21 to 10 s 21 . With these rates reduced to 0, the computations became unstable. A model was examined with the strain factor a for strain-dependent detachment (from states 1, 4, and 5) decreased from 2 pN 21 to a value of 0.5 pN 21 , making it close to the values obtained from measurements on myosin-actin binding [Nishizaka et al., 1995[Nishizaka et al., , 2000]. There were large decreases (30%) in association velocity and early dissociation velocity and a smaller increase in intermediate dissociation velocity. As a result, there was only a small decrease in dissociation velocity near the end of the P phase, and in some cases a linear fit was reasonable over the full travel of the dissociation transition. This change makes the results very different from the dissociation velocity results shown in Fig. 4 of Aoyama and Kamiya [2005].  Fig. 3, with addition of a linear "super-adhesion" elasticity of 0.12 pNÁnm 22 that restricts doublet surface separation beyond 34 nm. As the model contains one dynein per 12 nm length, this elasticity corresponds to 1.44 pNÁnm 21 per dynein, for comparison with values for FKY, usually less than 2 pNÁnm 21 per attached dynein. Qualitatively similar results were obtained when the separation distance at which super-adhesion elasticity was applied was varied from 30 to 38 nm, when the value of FKY varied from 0.2 to 4.0 pNÁnm 21 , or with the reduced stalk tilt models mentioned in the preceding section. Figure 7A covers one cycle of movement that shows bend initiation and propagation along the length of the doublet pair. Inclusion of interdoublet bending moment, as in Fig.  3A, is not essential for initial formation and propagation of a bend, but at least 30% of the full amount is required for regular cyclic behavior. The doublet pair model does not contain the elastase-sensitive shear resistance that constrains sliding in intact axonemes [Brokaw, 1980;Lindemann et al., 2005]. Consequently, with the full amount of interdoublet bending moment included in the model, a brief period of cyclic bending is usually terminated in a large bend that is limited only by the elastic bending resistance. This can be controlled by adding to the model a mathematical equivalent of an elastic shear resistance (not shown), but for the purposes of this article (Fig. 7) it has been controlled in a simpler manner, by reducing the interdoublet bending moment by a factor of 0.5. Figure 7B shows plots of doublet separation vs. length corresponding to each of the images in Fig. 7A. These plots show that the "super-adhesion" elasticity limits doublet separation to about 36 to 37 nm, which is more than sufficient for complete dynein detachment and inhibition of dynein force generation (Fig. 2B). Both separation at the distal edge of a separated region and reassociation at the proximal edge of a separated region are abrupt events that propagate together along the doublet pair as part of a propagated bend. The bends in separated regions are more prominent, with larger curvatures, than the bends in the opposite direction, between separated regions. The added curvature that is imposed by the interdoublet bending moment probably contributes to this difference. The bends between separated regions are the regions where dyneins are producing active sliding, which is responsible for bend propagation. Bend propagation velocity is about 300 mmÁs 21 .

Discussion Cyclic Splitting and Reassociation
The model used here for dynein shear force generation is similar to models used previously for simulation of flagellar movement. It is no surprise that it can also generate the shear forces needed to simulate the sliding that occurs in the experimental observations of cyclic splitting and reassociation by a doublet pair. This work now demonstrates that the same type of kinetic modeling produces dynein detachment and attachment that causes doublet separation and reassociation. The only novel feature is elastic compliance in the motor base, which forms its attachment to the A doublet. This elasticity increases strain on attached dyneins when the bending moments in the doublets are such as to produce doublet separation, as if pushed apart by a transverse force (t-force). This strain increases the rate constants for dynein detachment. The properties of the elastic attachment of the dyneins influence the results, but a range of properties is compatible with cyclic association and dissociation. Nonlinear elastic properties make the situation more realistic, by preventing dyneins from remaining attached over unrealistically long distances (more than 11 nm of stretch).
The dynein model used here is a simplistic model, which attempts only limited compatibility with information about the real structure of dynein motors, as reviewed in [Roberts et al., 2013]. It ignores knowledge that dynein in Chlamydomonas axonemes is a collective of many dynein varieties, possibly with functionally significant differences in distribution along the axoneme. One step toward reality has been taken by removing the assumption that a dynein attaches to a substrate site with a rigid, moment-bearing stalk. Instead, the dynein stalk can tilt, independently of its state, to produce the varied stalk tilts observed by Burgess [1995]. The assumptions that the power stroke is a linear motor producing movement exactly parallel to the A doublet, and that an adhesive force is derived from compliance that is exactly perpendicular to the A doublet, are idealizations that facilitate analysis. Deviations from these idealizations would probably be compatible with the experimental observations of Aoyama and Kamiya [2005] but have not been examined. The dynein model contains many parameters that must be specified, with limited guidance from chemical or structural data. Some combinations of parameters are better than others for replicating results available from the experiments. It would be possible to explore a much wider range of possible parameters, and gain more understanding of the operation of the doublet pair model. More data from doublet pair experiments would be needed to make this exploration worthwhile.
A more complete doublet pair model is introduced, which allows both doublets to bend in response to the forces generated by the dyneins. The models then show an overall bend comparable with the experimental example in Fig. 2 of Aoyama and Kamiya [2005], but the details do not exactly match the behavior of the experimental example. This is especially true when the bending moment resulting from interdoublet shear forces is included, as in Fig. 3A. In this example, near the basal end of the doublets when the doublets are not separated, there is an expected strong bending near the base. This results in a counterbend because of the viscous resistance to movement of the distal part of the doublet. The counterbend is difficult to see in

Bend Propagation and the "Geometric Clutch"
When both doublets are allowed to bend, it is possible to examine the relationship between bend propagation and doublet separation, as envisioned in the geometric clutch hypothesis for flagellar bend propagation [Lindemann, 1994a,b]. To do this, it is necessary to add to the model a replacement for structural connections that, in the absence of elastase digestion, normally constrain the separation of the doublets. This has been done by a simple mathematical construct that makes no attempt to represent real axonemal structures that can prevent excessive separation and at the same time allow interdoublet sliding in excess of 100 nm.
As shown in Fig. 7, this expanded model now demonstrates behavior anticipated by Aoyama and Kamiya [2005] on the basis of their occasional observation of doublet pairs with separations that propagated without much change in length of the separated region. The operation of the geometric clutch has usually been described in terms of dyneins being detached and deactivated when they are on the side of the axoneme where their sliding produced doublet separation, while doublets on the other side of the axoneme remain attached and activated. It is equally possible to describe its operation in terms of dyneins on the same side of the axoneme, in adjacent bends with opposite curvatures, remaining attached and activated. This view is necessary for interpreting a doublet pair, where there is no "other side." In either case, active sliding within a bend causes an increase in the magnitude of bend curvature at one end of the bend, and a decrease in curvature at the other end of the bend, which is equivalent to bend propagation. Propa-gation will be from base to tip in the majority of cases, where there is a strong shear resistance at the base. Tip to base propagation is predicted if the shear resistance at the tip is greater than the shear resistance at the base (Lindemann, 2007]. In the doublet pair model shown in Fig. 7, the bends in which active sliding is occurring are typically shallower and less obvious than the bends where the doublets have separated, but they are still bends. A new feature revealed by the modeling is the nonlinearity of doublet separation in response to curvature. This results in dynein force being turned on and off abruptly at the ends of a bend in which the doublets are separating. Previous modeling of the geometric clutch did not compute doublet separation, but just computed the transverse force (t-force) equivalent to the bending moments that are pushing the doublets apart and introduced a control relationship between t-force and dynein attachment probability. This modeling provided an initial verification of the geometric clutch hypothesis [Lindemann, 1994b[Lindemann, , 2002. It might now be improved by incorporation of the nonlinear switching revealed by the doublet pair computations. This is a useful interim approach that avoids the difficulty of computing forces and bending with variable interdoublet separation in a structure with more than two doublets.
The geometric clutch hypothesis received strong support from the experimental observations of Aoyama and Kamiya [2005], which indicated that control of dynein activity by doublet separation provides a sufficient explanation for oscillation in the two-doublet system. This interpretation was supported by modeling using simple mathematical equations for force generation and its control by doublet separation [Brokaw, 2009]. The present modeling takes the verification of the hypothesis a step farther. Realistic kinetic modeling of shear force and adhesive force generation, to produce and respond to doublet separation, continues to support the previous interpretation. Additionally, this modeling is used to demonstrate bend propagation when doublet separation is limited. The combination of experiments and theoretical modeling suggests that the geometric clutch may be a sufficient explanation for the role of bend curvature in controlling active sliding and flagellar bend propagation that was proposed many years ago [Brokaw, 1971]. Most importantly, the geometric clutch, first recognized by Lindemann [1994a,b], is not an added-on "control mechanism," but is an inherent property of a system that generates shear force between parallel doublets. The only requirements are a shear resistance that causes the shear force to generate bending, and adhesive forces that are not too strong to prevent separation and dynein detachment by the transverse force (t-force) resulting from shear force in a bent doublet. The dynein pair observations of Aoyama and Kamiya [2005], and their interpretation by computer modeling, demonstrate that these requirements are met. A different explanation is required for observations of oscillatory sliding in the absence of bending, such as Brokaw and Gibbons [1973] and Shingyoji et al. [1977]. Models that attempt to explain flagellar bend propagation entirely from an oscillatory sliding mechanism, such as Brokaw [1975] and Reidel-Kruse et al. [2007], could possibly be improved by allowing the geometric clutch to operate.

Methods
These computations of the movement of a doublet pair involve the repeated alternation of two steps.
Step 1 calculates the interdoublet forces produced by each dynein motor.
Step 2 uses these dynein motor forces to balance forces from elastic and viscous resistances and compute the bending of the doublets in a small time step, Dt. This information is then used to update the shapes of the doublets at the end of this time step.
Step 1 begins by using the current shape of each doublet to obtain the positions of the attachment site and potential attachment sites, of each motor, relative to the position of the motor. This information enables calculation of the strain on the elastic components of the motor, for each possible motor state, leading to calculation of the reaction rates for each transition between motor states. For each motor, the reaction rates are then used to calculate transition probabilities for transitions to each state accessible to that motor in the time interval Dt. A random number is generated, and then compared with the transition probabilities to determine the state of this motor at the end of Dt. The strains are then recalculated for this state of the motor, to obtain the motor forces to be used in Step 2. This stochastic procedure was developed [Brokaw, 1976[Brokaw, , 1995 and used for simulations of flagellar movement [Brokaw, 1999[Brokaw, , 2002. It is expanded here to incorporate a more complicated dynein model, with an adhesive force and a tilting stalk. These additions, and detailed definition of the dynein model shown in Fig. 1, are discussed below.
Step 2 uses methods developed previously [Brokaw, 1986b[Brokaw, , 2009] to solve the moment balance equation for the A doublet and obtain its movement as a function of time. As in that work, the dynein forces create longitudinal force in the A doublet that produces bending moment in the A doublet when it is bent. The modeling has been expanded for simultaneous solution of the moment balance equation for the B doublet. The first modification is to use the forces acting on the B doublet to calculate bending moments in the B doublet, just as is done for the A doublet. In addition, when both doublets can bend, the direct effect of bending moment that would normally bend both doublets together must be included. Half of this moment is included in the individual moment balance equation for each doublet. The bending moment is calculated as shear force between the doublets, multiplied by the normal doublet center-to-center distance of 60 nm. The use of this distance for the moment arm is an upper bound, which assumes that bending of a doublet is accommodated by compression and extension of its protofilaments, without shear between protofilaments. As significant shear force is only generated when the doublet separation is close to normal, neglecting any increased moment arm at increased separations is a reasonable approximation.
The doublets are constrained to bend in one plane, and no translation, rotation, or interdoublet shear is allowed at the basal end of the doublet pair. For numerical solution of the moment balance equation, the doublets are considered to consist of 112 segments with length 96 nm, and the integration proceeds with time steps of 0.002 ms. The total length of 10.752 mm includes 111 active segments and 792 dyneins (8 3 111 2 10% removed from tip region). One segment at the base allows no sliding and contains no dyneins.

Dynein Forces From Dynein Kinetics
Dynein kinetics are modeled using a 5-state mechanochemical ATPase cycle shown in Fig. 1 and augmented with more structural detail in Fig. 3 of Roberts et al. [2013]. This dynein model is similar to one previously used for modeling flagellar motility [Brokaw, 1999]. There are two strongly attached states, 5 and 1, before and after ADP release. Two detached states, 6 and 7, are added for "mechanical detachment" [Cooke et al. 1994] from the strongly attached states, bypassing the normal detachment from state 1 to state 2 resulting from ATP binding. Computations with the model actually use a total of 13 states. Another six states allow consideration of attachment to adjacent sites, one on each side of the site of current or most likely attachment, for each of the three attached states, 1, 4, and 5. Each dynein is evaluated at 0.002 ms intervals through time, to determine its current state and the force that it is exerting between the A and B doublets. Detailed specifications needed for dynein modeling are provided in Table I.
In its earliest forms, a dynein model was used to calculate shear force in the x direction, parallel to the A doublet surface. The power stroke was assumed to occur in the 2x direction and transfer its movement to the B doublet surface as if attached by a rigid, moment-bearing stalk that remained perpendicular to the A doublet. Electron-microscopic images of dyneins revealed variable stalk tilt [Goodenough and Heuser, 1982;Burgess, 1995]. Structural information about the dynein stalk [Ueno et al., 2008] suggests that it does not function as the "swinging cross-bridge" of a myosin motor. The Brokaw [1999] dynein model regarded the stalk as an extensible link, which also provided the elastic compliance of the motor in the x direction. This interpretation resulted in a nonlinear x-compliance Newer information about dynein stalk structure and its function in transmitting information [Gibbons et al., 2005] appears to be inconsistent with an extensible stalk, but does not fully explain how the stalk transmits force to the substrate doublet. In the model used here, the stalk is considered to have a fixed length (10 nm) and to be able to pivot at its attachment points to the motor (near the junction with the "buttress") and to a site on the B doublet, through a range of up to 6 0.45 rad from the y direction, perpendicular to the A doublet. The stalk tilt angle is not fully determined by the positions of the motor base and the B doublet site, measured parallel to the A doublet, because there are three variables: the stalk tilt, the strain in the y compliance of the motor, and the strain in the x compliance of the motor. The tilt angle is assumed to be the value that will minimize the total strain energy in these two compliances, and a simple one-dimensional optimization procedure is used to find the tilt angle for each calculation of the x and y strains and forces.
The elasticity of the y compliance needs to be nonlinear. A low resistance is needed for small separations, to allow buckling to initiate separation. A higher resistance is needed at larger separations, to prevent unrealistic extension of the base of an attached dynein motor. The compression When a motor is strained by forces in both the x and y directions, the force F used for calculating strain-dependent detachment rates is (F x 2 1 F y 2 ) 1/2 . c The ATP-dependent rate constants r 12 and r 72 for x 0 are calculated from r 21 and the energy level difference between unstrained state 1 and the next state 2, which depends on ATP concentration. r 12 and r 21 are reduced for x > 0 as with r 51 . d For motor base length y < 15 nm (compression): F y 5 2(y 2 15); For motor base length 15 nm < y < 16.5 nm: F y 5 FKY (y 2 15); For motor base length y >16.5 nm: F y 5 FKY (y 2 15) 1 2.5 FKY (y 2 16.5) 2 . resistance needs to be large enough to prevent buckling in the "wrong" direction.
Details are given in Table I.

Interactions Between Dynein Mechanics and Kinetics
Dynein mechanics will influence dynein kinetics in two distinct ways: 1. Beginning with the earliest model [Huxley, 1957], the rate of motor (cross-bridge) detachment was increased at the end of the power stroke, when the strain in the attached motor comes to 0 and reverses. This has been interpreted as a strain-controlled opening of the nucleotide binding site, allowing ADP to leave and ATP to bind [Smith and Geeves, 1995]. After ADP unbinding, detachment rate from state 1 to state 2 is determined by the rate of ATP binding, and not directly by strain. 2. An attached state, where the B end of the stalk is attached to a B doublet site, can be strained. The strain energy influences the attachment-detachment equilibrium; that is, the ratio between attachment and detachment rates. Energy levels that determine this ratio in the absence of strain must be specified for each of the attached states [Hill, 1974]. After one of these rates is specified, the other rate is determined by the ratio. Experimental work on strain-dependent breakage of myosin-actin rigor bonds [Nishizaka et al., 1995[Nishizaka et al., , 2000 suggested that the rate of strain-dependent detachment was proportional to e aF , where F is the applied force and a is an empirical constant, which was about 0.5 pN 21 for the myosin-actin measurements. The present modeling uses a similar formulation, r e aF , to calculate strain-dependent detachment rates for state 4 to state 3, state 5 to state 6, and state 1 to state 7, where r is the rate constant in the absence of strain. When strain is present in both the x and y directions, F 5 (F x 2 1 F y 2 ) 1/2 . For the weakly bound state 4, r has a large value (8000 s 21 ); for the strongly bound states 5 and 1, r has a small value (2 s 21 ). Even this value for the strongly bound state is higher than values for strain-free breakage of myosin-actin rigor bonds, reported to be less than 0.1 s 21 .

The Doublet Pair System
A value of 0.5 3 10 7 pNÁnm 22 has been used for the elastic bending resistance of a doublet. The value of 1 3 10 8 pNÁnm 22 determined for the elastic bending resistance of a relaxed sea urchin axoneme, after correcting for shear elasticity [Pelle et al., 2009], would suggest a value of 1 3 10 7 pNÁnm 22 for one of the nine doublets. The relatively high bend curvature measured during the normal bending cycle of a Chlamydomonas flagellum [Brokaw and Luck, 1983] suggests the use of a lower bending resistance than the value obtained for sea urchin flagella. The value of 0.5 3 10 7 pNÁnm 22 gives the doublets a realistic total bend when the model assumes eight dynein motors in each 96 nm segment, with average force of 2.7 pN/motor at 0 sliding velocity. This value would need to be adjusted to match different assumptions.
Interaction between the elastic bending resistance and the viscous drag on a doublet influences the rate of relaxation of doublets to their equilibrium conformations after separating completely from each other. A value of 2.16 3 10 29 pNÁsÁnm 22 has been used for the tangential viscous drag coefficient in simulations of experiments with sea urchin sperm flagella, such as Brokaw [1999]. Calculations based on the smaller diameter of a doublet and the different experimental conditions used with Chlamydomonas would suggest a value 0.7 times the value used for sea urchin sperm flagella. However, such calculations are only very rough approximations that ignore interactions between the doublets and other nearby surfaces. A factor of 0.58 instead of 0.7 has been used because it gives relaxation behavior leading to an association transition velocity close to the experimental measurements. When two doublets are close together, their total viscous resistance will be close to that of one doublet. Values used for the doublet drag coefficients are reduced by 0.55 when the doublet separation is less than 100 nm, but his adjustment has no noticeable effect. A drag coefficient ratio of 1.8 has been used, as in flagellar simulations.
A pair of partial differential equations, one for each doublet, must be solved simultaneously. For numerical integration of these moment balance equations, the length of the doublet is divided into N segments of length Ds. The doublet can bend at each of the joints, 1 to N 2 1, between segments. The distance of each joint from the base of the doublet is then measured by s. At any time t, the bend at each joint is represented by a curvature, j(s,t), which is the angle change at the joint divided by Ds. The array of curvatures at each joint completely defines the shape of the doublet at a particular time. The rates of bending at each joint, dj/dt, are the unknown quantities obtained by solving the moment balance equations. They are then integrated for one time step Dt to obtain values of j at t1Dt. A simple explicit integration, j(t 1 Dt) 5 j(t) 1 dj/dt Dt is unstable in the presence of large elastic resistances (a "stiff" equation), if dj/ dt is evaluated using moments at time t. Stability can be achieved by using a semi-implicit method, in which moments from elastic resistances are evaluated at time t 1 Dt and moments from external viscous resistances are evaluated at time t [Brokaw, 1985]. For a simple elasticity such as the bending resistance, EB, this is done by M(t 1 Dt) 5 EB j(t) 1 EB dj/dt Dt. In effect, this is including a viscous bending resistance equal to EBDt. Using external viscous resistances at time t, rather than at time t 1 Dt, introduces a slight error, but no instability. Moments resulting from shear elasticities, including the dynein-generated moments generated by strained elastic elements in the dynein motors, can also be treated implicitly. This is the most challenging part of the programming, because each Dt term involves a velocity between a point on the A doublet and a point on the B doublet, which depends on multiple values of dj/dt for both the A and B doublets.
For some computations, a super-adhesion elasticity was added, to produce an additional adhesive force, independent of dynein attachment, when doublet surface separation was beyond the distance required to detach the dyneins. This was done with a purely mathematical construct, by measuring distance to the B doublet surface, perpendicular to the A doublet, and applying an adhesive force along this line when the separation exceeded a threshold, such as 34 nm. A linear elasticity was sufficient. A similar compression resistance was used if the separation fell below 24 nm.