Introduction
Flagellar dyneins are a type of molecular motor located in the axoneme (see Fig. 1), where they drive the propagating bending motions of cilia and flagella, such as the beating of sperm tails (Roberts et al., Reference Roberts, Kon, Knight, Sutoh and Burgess2013). The canonical structure of motile cilia consists of nine pairs of doublet microtubules arranged radially around an additional central pair. The mechanical beating motion is generated by the coordinated activity of different isoforms of dynein motors arranged in linear arrays along the axoneme. The resultant effect of this coordination is to force adjacent pairs of microtubule doublet tracks to slide past one another (Summers and Gibbons, Reference Summers and Gibbons1971; Sakakibara et al., Reference Sakakibara, Kojima, Sakai, Katayama and Oiwa1999), causing the entire axoneme structure to bend (Riedel-Kruse and Hilfinger, Reference Riedel-Kruse and Hilfinger2007; Woolley, Reference Woolley2010). The tip of the dynein tail (see Fig. 1a) remains permanently anchored to a ‘cargo’ microtubule (see Fig. 1b), while the microtubule-binding domain (MTBD) (at the end of the stalk) repeatedly attaches to, exerts force on and detaches from the adjacent microtubule in a cyclical process, that leads to well-defined sliding motion between the two microtubule doublets. The directed force is applied as the motor transitions between two distinct conformational states, a transition known as the ‘powerstroke’, and therefore, the conformational states are known as the pre-powerstroke and post-powerstroke states.
Fig. 1. (a) An FFEA model of the pre-powerstroke (ADP•Vi, in red) and post-powerstroke (apo, in cyan) conformations of dynein-c. The tail domains of both conformations are superposed, showing the change in structure of the motor domain following the powerstroke. Microtubule-binding domain denoted MTBD. (b) An FFEA model of ADP•Vi dynein (in red) docked within the crowded environment of the axoneme (in white). The axoneme structure was obtained via cryo-ET. Microtubule denoted MT.
The kinetic cycle is driven by the binding and hydrolysis of ATP within the primary active site of the dynein motor, located within the head domain (Kikkawa, Reference Kikkawa2013). The pre-powerstroke has been observed in the presence of ADP-vanadate (ADP•Vi), which mimics the transition state of with ATP hydrolysis (ADP•Pi). In the post-powerstroke state the ligand is absent, so this is referred to as the apo state (Schmidt, Reference Schmidt2015).
Experimental studies of flagellar dyneins within the axoneme itself are difficult to perform due to the inaccessibility of the crowded environment. Nevertheless, the structure and flexibility of individual dynein-c motors have been determined, in the pre-powerstroke and post-powerstroke states, via cryo-electron microscopy (cryo-EM) (Roberts et al., Reference Roberts, Malkova, Walker, Sakakibara, Numata, Kon, Ohkura, Edwards, Knight, Sutoh, Oiwa and Burgess2012) and negative stain EM (Burgess et al., Reference Burgess, Walker, Sakakibara, Knight and Oiwa2003), respectively. Moreover, the 3D arrangement of these dynein motors and their tracks has been observed at ~10 Å resolution with cryo-electron tomography (cryoET) (Bui et al., Reference Bui, Sakakibara, Movassagh, Oiwa and Ishikawa2009). We have translated this multiscale experimental information into a continuum mechanical computational model, through the use of Fluctuating Finite Element Analysis (FFEA) (shown in Fig. 1), and used this representation to determine how the dynamics of an individual dynein motor are modified by its interaction with the complex environment of the axoneme architecture.
Methods
Fluctuating Finite Element Analysis
FFEA (Solernou et al., Reference Solernou, Hanson, Richardson, Welch, Read, Harlen and Harris2018) represents globular macromolecules such as proteins as continuum objects (Oliver et al., Reference Oliver, Richardson, Hanson, Kendrick, Read, Harlen, Sarah and Náray-Szabó2014), and has been previously used to model rotary ATPases, the ZipA/FtsZ membrane protein complexes, and cytoplasmic dynein (Richardson et al., Reference Richardson, Papachristos, Read, Harlen, Harrison, Paci, Muench and Harris2014; Lee et al., Reference Lee, Collins, Lin, Jamshad, Broughton, Harris, Hanson, Tognoloni, Parslow, Terry, Rodger, Smith, Edler, Ford, Roper and Dafforn2019; Hanson et al., Reference Hanson, Iida, Read, Harlen, Kurisu, Nakamura and Harris2020). Within the FFEA framework, proteins are approximated as viscoelastic solids (Wang and Zocchi, Reference Wang and Zocchi2011) which change shape and explore conformational space under the effect of thermal fluctuations. The 3D shape of the protein is represented using a volumetric finite element mesh, which can be constructed from cryo-EM or cryo-ET maps. Through the finite element method, the continuum equations of motion are transformed into a set of dependent linear algebraic equations for the set of nodes which define the mesh structure. Numerical integration of these equations then generates a trajectory describing how the shape of the protein changes with thermal fluctuations. The magnitude of these fluctuations at a given temperature is governed by the local material properties of the finite element mesh (within each tetrahedron), in particular Young's modulus, which is a required input parameter to the FFEA calculations.
The core continuum equation of motion defining an FFEA simulation is the Cauchy momentum equation. However, for the simulations performed in this work, it is appropriate to assume that the systems are overdamped, in which case our equation of motion reduces to:
where σ = σ elastic + σ viscous + σ thermal is the local stress tensor, composed of an elastic, a viscous and a thermal component. This summation of the viscoelastic components of stress constitutes a Kelvin–Voigt constitutive model, meaning that an FFEA object fluctuates about a predefined structure with no permanent deformation. Further details on the theory, the form of the stress tensor and the computational implementation of the method can be found along with the software release (Solernou et al., Reference Solernou, Hanson, Richardson, Welch, Read, Harlen and Harris2018).
Model construction
Cryo-EM structures of isolated dynein-c motors were used to construct a continuous 3D tetrahedral mesh of the motor, in both the pre-powerstroke and post-powerstroke configurations (Roberts et al., Reference Roberts, Malkova, Walker, Sakakibara, Numata, Kon, Ohkura, Edwards, Knight, Sutoh, Oiwa and Burgess2012) as shown in Fig. 1a. To correctly parameterise the dynein motors for subsequent calculations, we performed initial FFEA calculations of the isolated dynein motors (in the absence of the axoneme). We systematically varied Young's modulus of the stalk and the tail domains whilst keeping the Poisson's Ratio constant (Oliver, Reference Oliver2013), and calculated the different 3D conformational spaces explored by the motor as a function of the changing material properties. We compared the calculated range of angles and distances between the stalk and tail with those obtained experimentally (Burgess et al., Reference Burgess, Walker, Sakakibara, Knight and Oiwa2003), and selected values of Young's moduli that gave the best agreement with the experimentally determined range of motion, as shown in Fig. 2. Discrepancies in the positions of the peaks of the histograms, between the experimental and calculated distributions, arise from subtle differences between the equilibrium structures shown in the 2D negative stain EM image data, and those calculated from the 3D structures obtained via cryo-EM, i.e. the difference arises from the equilibrium structures provided by the two experimental methods.
Fig. 2. Angle and length distributions for dynein calculated from 2 μs FFEA simulations with Young's moduli that give the best fit to the experimentally determined ranges of angle and distance between the tail and stalk (Burgess et al., Reference Burgess, Walker, Sakakibara, Knight and Oiwa2003).
To accurately reproduce the experimentally observed dynamical behaviour, different Young's moduli were required for the stalk and the tail domains. This difference shows that for flagellar dynein molecules, the differences in internal atomistic structure between domains have a signature at the mesoscale that we can quantify within our continuum approximation. To prevent an abrupt change in the material parameters across the motor, we linearly interpolated the material parameters by distance across the motor domain, as can be seen in Fig. 3. The small aspect ratio of the motor domain should theoretically reduce any artefacts emerging from the interpolation. Additional input parameters used in the FFEA simulations, such as the internal viscosity of the proteins (Cellmer et al., Reference Cellmer, Henry, Hofrichter and Eaton2008) and the external viscosity of their surroundings can be found in Table S1 in the Supplementary Material.
Fig. 3. The Young's modulus distribution throughout each mesh for the ADP•Vi (Left) and Apo (Right) states. The tail and stalk regions are homogeneously assigned a modulus in accordance with the values specified in Table S3 in the Supplementary Material. The moduli of the motor-head region are determined by linearly interpolating between Young's moduli of the stalk and tail regions, to provide a smooth transition over the intermediate elements.
Conformational changes
We have seen that the constitutive model of stress within FFEA enables simulated proteins to fluctuate about a well-defined equilibrium position with no permanent deformation. This accounts for the equilibrium state of most naturally occurring proteins, which also have well-defined folded states about which they fluctuate. However, motor proteins such as dynein undergo radical conformational changes that are vital to their function. For dynein, the powerstroke is dependent upon ATP hydrolysis at the AAA1 region of the motor domain, which causes the linker region to undergo an approximate 90° rotation with respect to the rest of the molecule. In addition to the main results presented below, with FFEA we have the ability to perform such conformational changes directly whilst a simulation is in progress. We briefly demonstrate that functionality here.
As described previously, we have finite element meshes (node positions and tetrahedral elements) corresponding to each of the pre-powerstroke and post-powerstroke states of the dynein molecule. These meshes each have very different equilibrium shapes, and correspondingly also have different numbers of elements and numbers of nodes. To simulate the powerstroke itself, i.e. axonemal dynein switching between the two states, at some point in the simulation we need to replace one mesh with the other whilst keeping the current shape of the simulated molecule approximately fixed. To do this, we construct a mapping from the set of node positions of the pre-powerstroke state to the post-powerstroke state (and vice versa), accomplished as follows. We first create an FFEA simulation containing both a pre-powerstroke and post-powerstroke mesh. We identify a series of ‘equivalent points’ between the two structures (which are points that represent the same set of atoms in the two structures). We join each of these pairs of ‘equivalent points’ between the two structures with strong springs of zero equilibrium length (using linear elastic restraints within the FFEA simulation environment).
Given a sufficient number of such springs, we then run an FFEA simulation without thermal noise or steric repulsion between the two structures. The elastic restraints cause the two structures to overlap in a manner that minimises elastic energy, thus giving an overlapped configuration that is a compromise between the deformation of each individual structure.
This overlapped configuration allows us to define structural maps, such that the positions of the set of nodes of one structure can be formed as a linear interpolation of the positions of the nodes in the other structure (with a relatively small number of nodes requiring extrapolation). As there are different numbers of nodes in each structure, these maps take the form of non-square matrix transformations for both forward and backward transitions (from pre-powerstroke to post-powerstroke and vice versa). Since these maps allow the nodes from one mesh to be replaced by the nodes from the other mesh, they allow us to ‘switch’ from one mesh to the other when required, e.g. as determined from chemical rate kinetics.
Due to large thermal deformations and complex topological differences in the motor region, a transition is not always guaranteed to be possible using this mapping technique. Some mappings, for example, result in inverted mesh elements. In this case, the transition is rejected until the additional thermal motion of the current conformation renders the desired conformational change possible. Nevertheless, we have found that these issues with conformational mapping occur somewhat independently of the magnitude of the deviation of the molecule from the equilibrium structure. Further details of this methodology can be found with the original derivation (Hanson, Reference Hanson2018).
Figure 4 shows snapshots of a simulation performed using this conformational mapping technique to move through the conformational cycle, showing that the transition to the pre-powerstroke state enables the monomer to reach its next binding site location. To demonstrate this functionality, we provide a movie of an illustrative simulation of the chemo-mechanical cycle of axonemal dynein as Supplementary Material. Three distinct kinetic states are represented within the simulation, as shown in Fig. 4. When the motor is in the pre-powerstroke structural conformation it is coloured red, and when in the post-powerstroke conformation it is coloured cyan. Within the pre-powerstroke conformation, we further define a state of weak affinity for the MTBDs, shown in dull red, and a strong (binding) affinity for the MTBDs, shown as a vibrant red.
Fig. 4. A representation of the chemo-mechanical cycle of axonemal dynein simulated with FFEA. Conformational changes proceed from left to right, and then back to the start again.
During the FFEA simulation, we set the probability that the mapping will be applied (i.e. that a kinetic transition occurs) in accordance with a pre-defined transition rate. To approximately mimic the rates of the underlying biochemical transitions associated with dynein force generation, we assume that all conformational transitions occur at the same rate apart from the transition from the tightly bound post-powerstroke state to the unbound pre-powerstroke state. This rate is set to be extremely fast relative to the other transitions as once ATP has bound to the motor and triggered the release from the microtubule, no additional chemical triggers are required to trigger the subsequent conformational changes.
Results
Axonemal crowding
To quantify the minimal effect of crowding on the conformational flexibility of flagellar dynein within the axoneme, we constructed a model of the axoneme based on cryo-ET maps (see Fig. 1b), in which an individual dynein interacts only with the neighbouring microtubule tracks, and the neighbouring dynein motors are not represented. The cryo-ET maps show the positions of the motor domains, but the stalk and MTBD is highly flexible, as we have seen, and so is averaged out in the experimental image analysis.
We used a root mean squared (RMS) fitting procedure (see Supplementary Material) to dock our FFEA dynein mesh into the position occupied by dynein-c in the lower resolution structure, and produced a finite element mesh representing the intact dynein-c (in either conformation), and a surface mesh describing the adjacent (track) microtubule. The interaction between the microtubule and the MTBD is treated with a surface-surface Lennard-Jones-based interaction potential, representing Van der Waals (VdW) interactions from which the repulsive term forbids interpenetration of the finite element meshes. Specific binding sites for the MTBD are located at 8.1 nm intervals longitudinally on the protofilaments making up the microtubule surface, which are themselves staggered by 0.9 nm radially around the microtubule cylinder (shown explicitly in Fig. S2). We also used the attractive VdW interactions to generate the binding site adhesion for the MTBD.
During the chemo-mechanical cycle of the dynein motor, the MTBD (at the end of the stalk domain) in the pre-powerstroke (ATP or ADP•Vi) state firstly searches the microtubule surface for one of these binding sites. Following a tight binding phase, the conformational change into the post-powerstroke state occurs, followed by phosphate release and finally the release of ADP (Roberts et al., Reference Roberts, Kon, Knight, Sutoh and Burgess2013). To mimic the two biological situations for the pre-powerstroke and post-powerstroke states in the FFEA simulations, the favourable interaction potential between the MTBD and the track was set to a range of energies of the order k BT, so that weak binding (in which multiple binding and unbinding events could be observed) and strong binding states were explored. Using the material parameters obtained from fitting to the 2D EM flexibility data, we performed FFEA simulations and calculated the range of conformations explored by the motor in both the pre- and post-powerstroke states for each interaction energy. We also performed reference simulations in which the microtubule was non-interacting with respect to the motor, to measure the range of motion of dynein in the absence of crowding by its track. Further details of the FFEA parameters used to simulate individual dyneins and the generation of 2D representations of the 3D data can be found in the Supplementary Material.
Figure 5a shows the total volume of space explored by both states of the motor for each interaction energy between the MTBD and its interaction sites. In both the pre- and post-powerstroke states, an appreciably larger volume of conformational space is accessible to dynein when the repulsive interaction with the microtubule is switched off, and the motor is able to explore the volume that would be excluded by the track. The ability of the pre-powerstroke state (searching state) to locate subsequent binding sites is controlled by the reach of the dynein molecule whilst in this specific conformational state, and whilst confined within the axoneme environment. Figure 5b compares the probability density of the position of the binding domain for the pre-powerstroke state at an attractive track interaction energy of 1015Jm−4 projected perpendicular and parallel to the long axis of the microtubule. At this surface-surface energy density, which corresponds to a total favourable interaction energy of approximately 1 k BT, the simulations estimate the reach of dynein within the axoneme to be ~27 nm, which corresponds to the distance between the centres of the probability densities in the pre- and post-powerstroke states. In the pre-powerstroke state, conformational volumes associated with occupation probabilities of 30 and 95% contain 4 and 18 binding sites, respectively. For this single flagellar dynein interacting with its track, multiple binding sites located around, as well as directly along the track are accessible, implying that the reach of the motor need not be the same for each individual step. The ability of dynein-c to explore off-axial binding sites on a microtubule may contribute to the motor's ability to apply torque to the doublet, as has been observed in vitro (Schmitz et al., Reference Schmitz, Holcomb-Wygle, Oberski and Lindemann2000; Kikushima and Kamiya, Reference Kikushima and Kamiya2008). Tables 1 and 2 show the explicit quantitative results corresponding to these simulations.
Fig. 5. Probability density maps of MTBD spatial search as calculated from simulation trajectories of ADP•Vi dynein at different interaction energies,E int. Two viewing angles are shown for each energy: perpendicular to the axoneme (left) and down the centre of the axoneme (right). The head spends 95% of its time within the blue surface, 60% within the green, 30% within the orange and 10% within the red. Red arrows indicate the location of the motor head. (a) E int = 0 (b) E int = 1 k BT (c) E int = 10 k BT. Higher positive energies indicate more attractive interactions.
Table 1. The interaction potentials between the MTBD and microtubule-binding sites and their effect on binding sites exploration probability distribution
Higher positive energies indicate more attractive interactions.
Table 2. The interaction potentials between the MTBD and microtubule-binding sites and their effect on the volume explored by the MTBD
Higher positive energies indicate more attractive interactions.
To investigate the importance of representing specific binding sites on the microtubule, simulations in which the microtubule doublet surfaces were patterned with specific areas of higher interaction potential (the binding sites) were compared with a subsequent set of simulations in which the entire microtubule surface was uniformly attractive (no specific binding sites). No significant differences were detected in the explored-site counts, suggesting that the exact location of the binding sites may not be important in this case, presumably due to the large reach of dynein compared to the site separation on the microtubule surface.
Conclusions
Our FFEA simulations show that motion of the dynein stalk is impaired by the crowded environment of the axoneme relative to when it is in isolation on a surface, as is the case in the negative stain EM images (Burgess et al., Reference Burgess, Walker, Sakakibara, Knight and Oiwa2003). This is true even with a minimal representation, with only single dynein represented in the axonemal environment. The presence of other biomacromolecules, and the specific positioning and anchoring of the dynein motor relative to the track, modifies the probability distribution for the position of the MTBD compared to simulations in which hard-sphere repulsion is not included, and dynein can simply pass through the other objects within the axoneme. The implication is that quantities such as the affinity of the MTBD of flagellar dynein for microtubules, which has been measured by studying the binding of the relevant molecular fragments outside the native biological environment (Kamiya, Reference Kamiya, Keiko and Amos2012), may be modified in the context of the macromolecular organisation of the axoneme. More generally, and particularly whenever biomolecular complexes are located within complex 3D functional modules, it may be the case that biomolecular thermodynamics has an additional component that operates at the mesoscale, and which modulates atomistic interactions by way of the precise mesoscopic spatial and temporal organisation, in an analogous manner to the MTBD affinity for microtubules. Together with appropriate parameterisations at the mesoscale (Hanson et al., Reference Hanson, Iida, Read, Harlen, Kurisu, Nakamura and Harris2020), exploration of this new layer of biological complexity by combining conventional biophysical measurements, such as molecular binding assays, with emerging data from experimental techniques, such as cryo-ET and super-resolution microscopy (Lučić et al., Reference Lučić, Rigort and Baumeister2013; Gray et al., Reference Gray, Harlen, Harris, Khalid, Leung, Lonsdale, Mulholland, Pearson, Read and Richardson2015), will be instrumental in furthering our understanding of the physical mechanisms underlying the operation of biomolecular machines across multiple length and time scales.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033583520000062.
Acknowledgements
We would like to thank Dr Takashi Ishikawa, who provided the density map for the two microtubule doublets built from the data in reference (Bui et al., 2009), Stan Burgess and Michelle Peckham for helpful discussions and Molly Gravett for reading the manuscript. RA Richardson received an EPSRC DTP studentship, and additional support was made by the EPSRC through a DTA to BS Hanson. Simulations used the ARC High-Performance Computing facilities at the University of Leeds, UK.
