A mathematical understanding of how cytoplasmic dynein walks on microtubules

Cytoplasmic dynein 1 (hereafter referred to simply as dynein) is a dimeric motor protein that walks and transports intracellular cargos towards the minus end of microtubules. In this article, we formulate, based on physical principles, a mechanical model to describe the stepping behaviour of cytoplasmic dynein walking on microtubules from the cell membrane towards the nucleus. Unlike previous studies on physical models of this nature, we base our formulation on the whole structure of dynein to include the temporal dynamics of the individual subunits such as the cargo (for example, an endosome, vesicle or bead), two rings of six ATPase domains associated with diverse cellular activities (AAA+ rings) and the microtubule-binding domains which allow dynein to bind to microtubules. This mathematical framework allows us to examine experimental observations on dynein across a wide range of different species, as well as being able to make predictions on the temporal behaviour of the individual components of dynein not currently experimentally measured. Furthermore, we extend the model framework to include backward stepping, variable step size and dwelling. The power of our model is in its predictive nature; first it reflects recent experimental observations that dynein walks on microtubules using a weakly coordinated stepping pattern with predominantly not passing steps. Second, the model predicts that interhead coordination in the ATP cycle of cytoplasmic dynein is important in order to obtain the alternating stepping patterns and long run lengths seen in experiments.

Cytoplasmic dynein 1 (hereafter referred to simply as dynein) is a dimeric motor protein that walks and transports intracellular cargos towards the minus end of microtubules. In this article, we formulate, based on physical principles, a mechanical model to describe the stepping behaviour of cytoplasmic dynein walking on microtubules from the cell membrane towards the nucleus. Unlike previous studies on physical models of this nature, we base our formulation on the whole structure of dynein to include the temporal dynamics of the individual subunits such as the cargo (for example, an endosome, vesicle or bead), two rings of six ATPase domains associated with diverse cellular activities (AAAþ rings) and the microtubule-binding domains which allow dynein to bind to microtubules. This mathematical framework allows us to examine experimental observations on dynein across a wide range of different species, as well as being able to make predictions on the temporal behaviour of the individual components of dynein not currently experimentally measured. Furthermore, we extend the model framework to include backward stepping, variable step size and dwelling. The power of our model is in its predictive nature; first it reflects recent experimental observations that dynein walks on microtubules using a weakly coordinated stepping pattern with predominantly not passing steps. Second, the model predicts that interhead coordination in the ATP cycle of cytoplasmic dynein is important in order to obtain the alternating stepping patterns and long run lengths seen in experiments.
integrating the temporal dynamics of the individual components that include the cargo, tail domain, AAAþ rings and MTBDs. For simplicity, we formulate our modelling on a one-dimensional microtubule, leaving extensions to multi-dimensions for future studies. Similarly, modelling of multiple dyneins walking on microtubules [7] is omitted and forms part of our future studies. In §5, stochasticity is introduced into the model to account for the random binding of ATP to either of the two motor domains, and numerical simulations for the model equations are presented. Within this section, qualitative and quantitative agreements with some experimental observations are discussed, and the effect of interhead coordination is explored. Backward stepping, variable step size and dwelling are further modelled, and numerical simulations exhibit this stepping behaviour. Furthermore, we make predictions amenable for experimental manipulations. Finally, in §6, we discuss the implications of our modelling to understanding mechanisms for dynein-mediated transport.

Experimental observations
Experimental studies (using total internal reflection fluorescence, X-ray crystallography and highresolution cryo-electron microscopy) on how cytoplasmic dynein motors move along microtubules transporting cargo to the nucleus can be subdivided into two parts. Those that focus on the single molecule motility properties of dynein [3,21,36,37] and those that focus on the mechanical structural dynamics of dynein [2,20,21,35,38]. Furthermore, experimental data both in vivo [39][40][41] and in vitro 3 [2,4,42] are performed and obtained on different species (e.g. yeast [3,35], Dictyostelium discoideum [2,4,7,20,23], human [38,42], Saccharomyces cerevisiae [21,36,37], etc.). Moreover, studies are carried out either in one-dimension [37] or two-dimensions [4,36,37]. In the work by Qiu et al. [37], twodimensional particle tracking shows dynein's two motor domains can step both alternatively and nonalternatively in time and either passing or not passing in space. One-dimensional tracking results are then extrapolated from the two-dimensional data through a projection operator in the direction of motion along the microtubules axis. Given these different experimental conditions, it is therefore a significant challenge to come up with a single mathematical model that can capture all these processes. Given that the single molecule motility dynamics are part of the whole dynein structure, we therefore propose to study the whole dynein structure with the rationale that different experimental conditions can either be modelled through parameter variations and functions or through appropriate extensions of the model to take into account other processes associated with dynein transport that are not the subject of our study. Having this in mind, we therefore present some of the most recent experimental observations of cytoplasmic dynein with an eye to making comparisons with the model where appropriate. Moreover, the structural mechanism of cytoplasmic dynein's processive stepping along MTs is unclear and sets the motivation for this study.
Recent experimental studies on the structure of yeast and Dictyostelium discoideum dyneins include works by Bhabha et al. [35], Carter et al. [18] and Schmidt et al. [21], which allowed for detailed visualization of the AAA domain and linker movements. On the other hand, studies on yeast and Dictyostelium discoideum dyneins show that the replacement of the tail with glutathione S-transerase yields a simpler dimer that still processively steps along MT [7]. It is known that perhaps the most striking feature of stepping dynein is the huge flexibility between the ATPase domain and the trackbinding domain, which is in contrast to kinesin and myosin motors [7]. In order to provide quantitative data on dynein's processive stepping along MTs, Dewitt et al. [36] and Qiu et al. [37] tracked fluorescent-tagged yeast cytoplasmic dynein in two-dimensions, and their studies suggested uncoordinated stepping pattern by the two heads but that they also must communicate (i.e. coordinate) as the properties of dimerization to MTs are different from those of monomers [7]. Statistically, their studies found that 74% of dynein steps were taken by each of the two heads alternating in time (i.e. the motor domains' relative temporal behaviour) and that 83% did not pass each other (in terms of their relative spatial behaviour), suggesting that dynein may move predominantly by passing rather than in an alternating fashion [37]. This is in contrast to the processive kinesin and myosins which walk handover-hand [4,[43][44][45][46]. Both papers also found that the leading head was more likely to be to the right of the lagging head along the direction of movement [36,37]. Furthermore, experimental observations show that dynein has a variable step size, with the majority of steps being 8 nm in distance [36,37,[47][48][49]. It must be noted, however, that this reflects the position of the tail of dynein rather than the motor domain, and further investigations have shown that the motors move with a usual step size of around 16 nm [36,37,48]. Moreover, we note that dynein steps are not always parallel to the microtubule and usually have off-axis components [10,36,37,48], and dynein can also take backward steps [47,48]. Carter et al. [50] proposed that the stalk acts as a tether in the stepping process and that the MTBD determines the direction of the step, while Redwine et al. [51] proposed that conformational changes in the MTBD lead to movement in the linker domain and hence displacement of the MTBD.
We note that the experimental studies by DeWitt et al. [36] and Qiu et al. [37] revealed the stepping behaviour of yeast cytoplasmic dynein; however, throughout this paper, we will also consider mammalian cytoplasmic dynein; hence, the model is not restricted to specific species. Furthermore, due to the fact that at the stalk-stalkhead junction, the hinge is located close to the MT surface (twodimensional geometry), the dynein head swings over a wide range of approximately 20 nm compared with approximately 8 nm spacing between the binding sites on the MT. Studies by Imai et al. [7] suggest that experiments such as those by DeWitt et al. and Qiu et al., where fluorescent tags are attached to the heads for stepping studies, may not reliably report the position of the MT-bound stalkheads. During the processive stepping, the flexibility between the ATPase and track-binding domains may allow for the stalkhead to detach from its partner motor with greater freedom to explore the MT surface for locating its next binding site. Hence, these studies provide a structural basis for a wide range of step sizes (variable) seen in dynein stepping studies [7]. As a proof of concept, we will nevertheless compare our results to those of DeWitt et al. and Qiu et al. (with the caveat above and noting also that our model is formulated in one dimension) as our theoretical study is a first stepping stone in modelling the integrated dynein structure. The framework can easily be applied to specific dynein stepping studies when quantitative experimental data are available. For example, by including a linker into the modelling, results could be compared quantitatively with those obtained in studies by Cleary et al. [4]. While these experimental studies enable detailed understanding of dynein's structure [2,4,7,18 -21,35] and transport mechanism [3,34,36,37,42,47,48,52]; to our knowledge, few models have been developed to describe such observations [53][54][55]. Our results bridge this gap, by presenting a robust integrative mechanical and stochastic model describing the stepping behaviour of cytoplasmic dynein.

Overview of current mathematical models describing transport processes for cytoplasmic dynein
Recently, there has been a notable increase in mathematical models studying endocytosis and cytoplasmic dynein and these include models proposed by Ashwin et al. [56], Smith & Simmons [57], Šarlah & Vilfan [54], Mukherji [58] and Tsygankov et al. [55,59]. For example, the Smith and Simmons model [57] allows for motion of dynein along the microtubules when cell organelles and vesicles, referred to as particles, are attached and freely diffuse when they are not. Under this framework, they consider particle densities in one dimension, described by reaction-diffusion-transport equations. This model helps us to understand the macroscopic behaviour of endosomes and not the particular mechanisms of dynein. In the model proposed by Ashwin et al. [56], they consider a single microtubule for which the motor protein dynein moves to the minus end (i.e. towards the nucleus) when bound and is carried by the motor protein kinesin to the plus end (towards the cell periphery). They assume that right and left moving motors pass without interaction, but there is an exclusion principle enforcing that a motor can only move forward if the site ahead is free of motors of the same type. They discretize the microtubule into two tracks and use a mean field approximation and further simplifications. This model describes the behaviour of a population of dynein; it does not consider a more detailed model involving a single dynein within the transport process and therefore is not able to quantify the temporal behaviour of the individual components of the dynein structure. Single dynein models have been considered by Mukherji [58], Tsygankov et al. [55,59] and Šarlah & Vilfan [54]. Mukherji [58] and Tsygankov et al. [59] study the mechanochemical cycle of dynein which is essential for understanding dynein's behaviour. An extension to the model by Tsygankov looks at the bending energies of dynein using Langevin equations and couples this to the biochemical reactions modelled previously [55,59]. Šarlah and Vilfan propose a winch model for cytoplasmic dynein which couples an elastomechanical model to a kinetic model of the ATPase cycle. For the elastomechanical model, they consider elastic energies within the complex, interaction between the two motors and work done against external load. Monte Carlo methods are then used to find the shapes of the complex with minimum energy. In this work, we propose an alternative framework; our aim is to derive a model that studies dynein's progress along the microtubule over time as opposed to mean run lengths and velocities, taking a mechanical approach with the long-term aim of modelling the mechanical effects of mutations on dynein. A similar approach has been studied for kinesin by Hendricks et al. [60] and by us in a previous work by Crossley et al. [53]. We will study the whole structure, looking at the positions of the cargo carried by dynein, the tail domain, AAAþ rings and MTBDs comparing our results to data from different experiments tracking single components of the transport process. We will also consider dynein in general, allowing the model to be applied later to particular dynein species through the use of parameter variations. Unlike this current work, the previous study was devoid of any statistical analysis which forms the bulk of the current modelling approach. The main contributions of this study are the stochastic multiscale modelling, as opposed to the use of continuous functions to model binding and ATP force previously studied, and the introduction of the tail component. Furthermore, we model for the first time variable stepping, including heads being able to move independently, not in a strictly coordinated pattern. Other processes such as variable steps, backward stepping and dwelling times are also modelled for the first time.  The coordinates x A and x D represent one head domain of dynein with the coordinates x B and x E representing the other head (figure 2). We model the microtubule as a one-dimensional line with binding sites 8 nm apart; we only consider motion along this line. We make the following assumptions:

Derivation of the mechanical model
-The mass of the cargo remains constant and is modelled as a sphere with small Reynolds number. This is a significant assumption for experiments in vivo; however, it is applicable to in vitro experiments with beads. -Any regulators of cargo binding, such as dynactin, are modelled as part of the cargo. -The tail domain is modelled as two identical springs, from the AAAþ rings, connected to a sphere with small Reynolds number and constant mass. The linker is modelled as part of these springs. The binding between the tail and the cargo is modelled via another spring connecting the tail domain to the cargo. -The AAAþ rings are identical and modelled as spheres with small Reynolds number whose masses remain constant. -The stalks are modelled as two identical springs. We model the strut or buttress as part of this spring. -The MTBDs are identical and modelled as spheres with small Reynolds number whose masses remain constant.
See figure 2 for a schematic diagram illustrating the whole structure on which the mathematical model is based and table 1 for a list of parameter values. We make four simplifying assumptions that will be relaxed in future studies (see remark 4.1): -The spring between the cargo and tail domain is parallel to the microtubule. Remark 4.1. The simplifying assumption of fixed angles means that the AAAþ rings and cargo will move according to the extension and relaxation of the springs horizontally. This is an appropriate assumption for the model while we remain in one space dimension but will need to be considered when moving to higher dimensions. It is likely that there is some rigidity within the complex with regard to these angles, with the main variation arising from the conformational change under ATP MTBDs drag Figure 2. A schematic diagram of the mechanical model (adapted from [53] It must be noted that structural studies show that both stalks are tilted towards the plus-end and more or less in parallel orientation to each other (rather than pointing towards each other as depicted Table 1. Dimensional parameters and the primary values used in the mathematical model. The drag coefficients are given by g i ¼ (6phR i ) MDa/ns for i ¼ C, T, M, S with h and R i given below. The binding sites are described by p kþ1 ¼ ( p k þ 8) nm where p 0 ¼ (L C þ L T 2 4) nm, with L C , L T given below. in figure 2). The mathematical framework can be modified easily to take into account this particular structure with no changes in the numerical results and model predictions (see the electronic supplementary material for further details on the model reformulation and the corresponding numerical results). Hence, our studies confirm similar results if stalks are assumed to be tilted towards the plus-end with angles ranging between 41.98 + 13.78 [2,7,29,34,64,65]. In order to account for the angles of the dynein's off-axis steps, it is necessary to consider a two-dimensional model. The two-dimensional model will allow us to investigate whether dynein has a preferential stepping behaviour, either right or left on the microtubule surface. Such an analysis is not possible within the one-dimensional set-up proposed in this study. To proceed, using Newton's Second Law we study the net forces acting on the system. For the cargo, there is a spring force, viscous drag and an external force acting on it. By Hooke's Law, we take the spring force to be: where K C is the spring constant and L C is the natural length. We obtain the viscous drag by Stokes' Law: where the damping coefficient g C ¼ 6phR C with h the viscosity and R C the radius of the cargo. For completeness, we include an external force F C that is exerted on the cargo, although throughout the model this is assumed to equal zero. Therefore, the equation of motion for the cargo can be modelled by The equations of motion for the tail domain and AAAþ rings can be derived similarly. Therefore, we obtain the following system of ordinary differential equations for the cargo, tail and AAAþ rings, respectively: We wish to model the mechanics of ATP hydrolysis on the motor domain of dynein. The binding of ATP occurs randomly and is followed by microtubule release of the corresponding MTBD and a recovery stroke towards the next binding site [16,24]. Hence, we will assume that there are two MTBD states: -Bound: This is defined to be when the MTBD is bound to the microtubule and hence is stationary.
-Unbound: Defined to be when the MTBD is unbound from the microtubule and undergoing the recovery stroke towards the next binding site.
It is hypothesized that ATP hydrolysis induces a conformational change in dynein, potentially causing a 378 kink in the stalk [19]. Hence, for the unbound state the conformational change is modelled by a dashpot and spring acting solely on the MTBD (figure 3) [66]. It is assumed that this force is independent of the particular interval on the microtubule, defined by x [ [ p k , p kþ2 ], and is identical for the two head domains. Binding sites are taken to be p 2k for MTBD D and p 2kþ1 for MTBD E with k ¼ 0, 1, 2, . . . and p 2kþ1 2 p 2k ¼ 8 nm, binding sites are 8 nm apart on the microtubule with each MTBD binding to distinct binding sites that are 16 nm apart. The current model is one dimensional and hence it is assumed that this force acts only in the horizontal direction. The force rsos.royalsocietypublishing.org R. Soc. open sci. 5: 171568 produced by the dashpot is proportional to the velocity and the spring force is proportional to the displacement, hence where g ATP and K ATP are parameters determining the size of the ATP force, with estimated values given in table 1. These parameters are estimated by trial and error. The parameter L ATP represents the unstressed length of the spring and is taken to be the step size of the head domain. Here, we use a fixed step size of 16 nm; however, in §5.4 we model variable step sizes in order to represent more faithfully experimental observations which show dynein stepping in a variable fashion. If MTBD D is in an unbound state and MTBD E is in a bound state, then the equations of motion can be shown to be given by is the binding site that MTBD D was bound to at time t ¼ t i . The equations are similar for when MTBD E is in the unbound state and MTBD D is in the bound state: 12Þ Here, we are assuming some inherent coordination between the two MTBDs to keep the motor attached to the microtubule as one motor is unable to bind ATP, while the other is detached. The MTBDs are assumed to become unbound once the corresponding AAAþ ring binds ATP. This occurs randomly and the transition between states is explained below. The model is extended to include dwelling between steps, backward stepping and a variable step size in §5.4.
Remark 4.2. Fixing the step size to 16 nm with predetermined binding sites is a strong assumption on where the MTBDs can bind. MTBDs are restricted to binding to specific binding sites on the microtubule due to the position of tubulin and cannot bind to a site that another MTBD is already bound to. The displacement of the MTBD, under a conformational change during the ATP cycle, has been suggested to be close to the 16 nm step size [19], with this being the predominant step size in these studies. Other step sizes have been recorded alongside off-axis displacement. For simplicity, we will consider the simplest model of dynein stepping in one space dimension. Variable steps sizes are considered in §5.4, while two-dimensional stepping is left for future studies.

Continual stochastic stepping
To model the continual stepping by dynein over a microtubule, stochasticity is introduced to the model via the randomness in which an AAAþ ring binds ATP and hence an MTBD becomes unbound. We assign the values P D (P E ) to the probability that MTBD E steps given that MTBD D (E) stepped previously and the maximum separation distance d that can occur between the MTBDs is defined.
be a random vector where q i is from the uniform distribution on the interval (0, 1). If the maximum separation between the MTBDs has been exceeded, then it is assumed that the rearward head steps; else, given that MTBD j stepped previously, if q i , P j , then MTBD E is set to be in the unbound state (i.e. unbound from the microtubule and undergoing the recovery stroke) and MTBD D is set to be in the bound state (i.e. bound to the microtubule). Otherwise, we assume that the MTBD D is in the unbound state and MTBD E in the bound state. Hence, we can define a step function h E given by This does assume some form of coordination between the two head domains of dynein as only one head will step during each time interval, but it does not enforce coordination of the stepping pattern itself if the head domains are allowed to separate past consecutive binding sites. The rearward head always steps if the two head domains become too far apart. This assumption reflects the existence of a linker that plays a critical role in gating dynein stepping behaviour thereby modelling tension-dependency at high interhead separation during the processive stepping [4]. In future studies, it might be worth introducing a model specifically taking into account how the linker gates dynein stepping behaviour [4]. The system of ODEs is therefore given by  In this model, we only consider continual stepping; therefore, we fix the size of the time interval for each step, T Step ¼ t iþ1 2 t i , and hence T Final will depend on the time interval T Step and the number of steps N. Therefore, the stepping rate of the motors is predetermined. This assumption is relaxed in § §4.4 and 5.4 where independent and random dwell times are introduced to the model, respectively.
Remark 4.4. The binding sites are predetermined. The initial binding site p 0 is assigned a value and all binding sites are taken to be 8 nm away from the previous binding site. For each time step, the binding site is updated by taking the next binding site of the unbound MTBD. For example, if an MTBD is unbound on [t i , t iþ1 ] and bound to p k at time t ¼ t i , then the binding site will be updated to rsos.royalsocietypublishing.org R. Soc. open sci. 5: 171568 It must be observed that the time continuity of the model is reflected and embedded in the time-variables associated with the domains. These variables are monitored (once computed) if they are located at the discrete steps or not. Variable step sizes are explored in §5.4; however, they are restricted to multiples of 8 nm to ensure that they can only bind at a specified binding site on the microtubule.

Non-dimensionalization
To non-dimensionalize the model, let The non-dimensionalized coefficients of the acceleration terms turn out to be small and the dynamics are dominated by the viscous drag [53]. Hence, neglecting the small coefficients of the second derivatives, we obtain the following non-dimensional system: The non-dimensional parameters are given by

Initial conditions
We prescribe initial conditions as follows: MTBD D and E are taken to be at binding sites p 0 and p 1 ¼ ( p 0 þ 8) nm, respectively. The cargo is taken to be at the origin and the tail component is set to be at its natural length L C from the cargo. The AAAþ rings are taken to be at the same point midway between the MTBDs, at a distance of the natural length L T from the tail. Therefore, the initial conditions are set to be The non-dimensional initial conditions are given by

Independent stepping
It has been suggested in previous studies that interhead coordination is important to the stepping mechanism of two-headed cytoplasmic dynein [67]; therefore, we explore the significance of this coordination by considering the resultant behaviour if it is disrupted, i.e. if the two head domains step independently. The dwell time before the binding of ATP for each motor domain is modelled by the exponential distribution (see remark 4.5). We assume that these waiting times for each motor domain are independent of each other and are given by MTBD E with q i D and q i E taken from the exponential distribution with mean dwell time m. The system continues to be modelled by equations (4.10)-(4.13), with different stepping functions to h D , h E in equations (4.14) and (4.16). For MTBD D, we assume that it steps after q i D ns, hence we define the following step function: where t i and t iþ1 are the times when MTBD D binds back onto the microtubule after stepping with t 0 the initial time. The stepping function for MTBD E can be defined similarly: with t j and t jþ1 the times when MTBD E binds to the microtubule. Here, t j denotes different time intervals to t i . Therefore, the following model system of ODEs can be derived as follows: Remark 4.5. Experimental observations suggest that the dwell times of dynein can be approximated well by an exponential distribution with an average dwell time of 2 s, i.e. 2 Â 10 9 ns, per ATP cycle [48]. It is currently assumed that the dwell times are identical; however, differences in mean dwell times could be explored in future work.
We take a multiscale approach when non-dimensionalizing the model system, using one fast timescale for the stepping and one slow timescale for the dwelling. For the dwelling interval, we non-dimensionalize as above with t c ¼ m and for the stepping intervals we take t c ¼ m C /g C . For the sake of brevity, details of the non-dimensionalization are omitted here (see the electronic supplementary material for details).

Numerical experiments
The scheme is implemented in MATLAB for N stepping intervals of [0, T Final ] with non-dimensional end time T Final ¼ 10 8 and N ¼ 100 using the solver ode45 [68]. ode45 is one of several solvers for integrating a system of non-stiff ordinary differential equations given appropriate initial conditions and is based on Runge -Kutta time-integrators. For further detailed description and implementation, we refer the interested reader to consult MATLAB MathWorks [68]. The initial conditions are given by For the initial step, it is assumed that MTBD D is in an unbound state and MTBD E is in a bound state. For each following step, a random number q i is generated from the uniform distribution on the interval (0, 1) and the initial conditions are given by the values from the previous simulation: . From here onwards, we refer to trajectories, the physical loci or paths taken by each x i (t), and these represent the distances travelled in time. Also these could be referred to as positions of the components as a function of time.
Remark 5.1. In all our numerical simulations (unless stated otherwise), in the absence of explicit modelling of the tension generated by the linker to gate dynein stepping behaviour, we impose a maximum interhead separation distance of 48 nm.

Stochastic stepping with limited coordination
Initially, we assume that the motor domains bind ATP at random when they are both attached to the microtubule; therefore, we take P D ¼ P E ¼ 50%. This entails that the two head domains will not be highly coordinated in terms of their ATPase cycle, although they will experience some coordination, in terms of attachment to the microtubule, as we assume that only one motor domain can detach at a time. The results show a mixed stepping pattern for both the MTBDs and AAAþ rings with both not passing and passing stepping patterns present (figure 4d,e). This matches experimental observations (onedimensional projections of two-dimensional experiments) on yeast cytoplasmic dynein, labelled at the AAAþ rings [36,37]. Here, we are able to compute the trajectories of the AAAþ rings and MTBDs, which is not yet achievable in experiments, as tagging functional MTBDs is technically challenging. The tail domain also moves with a stepping profile, as seen in experiments on dynein labelled at the tail domain ( figure 4c,f). The cargo moves along the microtubule with increasing velocity, which becomes oscillatory at longer times once the dynein has settled into a stepping behaviour ( figure 4a,b). By computing the solutions over a larger interval, with end time t Final ¼ 10 9 and N ¼ 1000, the velocity of the cargo reaches a relative plateau where it stops increasing over time and oscillates within a small band ( figure 5a,b), matching observations by Garrett et al. [10]. The velocity of the cargo increases over time before reaching a plateau due to the fact that dynein starts to move from a stationary position and therefore must pick up speed. We do not impose an external force at the initial phase of the stepping process.
In order to make statistical comparisons with experimental observations, we compiled data from 1000 simulations of the model and then took averages over these different realizations. Observations by Qiu et al. [37] show that approximately 83% of steps did not pass each other and in our simulations we have an average of 84.71% steps not passing when d ¼ 48 nm (figure 6a, table 3 6b, table 3). This may be due to the randomness in the model where the probability of stepping is independent of which head stepped previously.  [37], use dimerized yeast dynein. We have therefore also looked at a reduced version of the model for a dimerized dynein motor with no cargo and we get similar results for the stepping pattern and trajectories (see the electronic supplementary material for details).

Extensive interhead coordination
If dynein uses a more extensive form of interhead coordination, the probability that each MTBD steps will depend on the previous step. Therefore, the impact of dependent stepping probabilities on the model is investigated by taking P D = P E . It is assumed that the probability that MTBD E steps increases if MTBD D stepped previously and decreases if MTBD E stepped previously. By taking P D ¼ 70% and P E ¼ 30%, the results show the same mode of stepping to previous results, with a mixed stepping pattern of 84.0% not passing steps reflecting experimental observations of 83% (figures 7 and 8a). However, in comparison to our previous results, these results also resemble

Independent stepping
We now relax the assumption that there is coordination between the head domains and that they step independently. Owing to the independence of the two MTBDs, both MTBDs could become detached from the microtubule, if this occurs then the simulation is terminated and the number of steps and the run length are recorded. Initially, we consider forward stepping with a fixed step size of 16 nm. We consider a maximum of N ¼ 100 steps with a mean dwell time of m ¼ 2 Â 10 9 ns for each head domain. Numerical simulations are run in MATLAB using the solver ode15s for the dwelling period and ode45 for the stepping intervals [68]; example profiles are given in figure 9, and the mean percentage of steps passing and alternating are given in figure 10. The MATLAB solver ode15s is employed here to provide numerical solutions to a system of stiff ordinary differential equations (and in practice as well as systems of differential-algebraic equations (DAEs)), unlike ode45 previously used [68]. By analysing the stepping behaviour of this model, we see that for larger values of K ATP , in the 400 -1000 MDa ns 22 range in table 5, we achieve 83.63% to 86.97% not passing steps on average, which is close to the 83% seen in experiments. However, all values of K ATP in table 5 give much lower values for the average percentage of alternating steps than those seen experimentally (predominantly around 49% compared to 74%). This suggests that this independent form of stepping cannot account for the alternating stepping patterns seen in experiments and hence there must be some form of coordination acting between the two head domains to account for this behaviour. Remark 5.3. We note that by allowing the head domains to step independently, they are able to diverge considerably implying larger interhead separation distances which might not be biologically realistic (e.g. figure 9c,d ). This shows the need to explicitly model the linker which has been shown to  gate interhead coordination at larger separation distances (and therefore able to bring the domains closer to each other) [4]. Alternatively, modelling internal forces that are known to influence the stepping behaviour may bring the complex back together [47]. Both of these remedies are the subject of our future studies.
Variations in the values of K ATP show that run lengths are highly dependent on this parameter. We see that for K ATP ¼ 10 MDa ns 22 and K ATP ¼ 100 MDa ns 22 the mean number of steps in a run is less than one, suggesting that predominantly the run is terminated before the first step can be completed. Processivity is therefore dependent on the value of K ATP . Although this parameter cannot be directly measured in experiments as it is an approximation of the effects of the ATP force, it suggests that if the ATP cycle of the detached head domain is not completed quickly enough, then an uncoordinated detachment of the attached MTBD is likely to occur and hence the run will be terminated after fewer steps. Taking K ATP ¼ 500 MDa ns 22 gives a mean number of steps of 33.50 and mean run lengths of 275.95 nm for the cargo and 276.69 nm for the tail domain (tables 5 and 6). Although this gives the highest run lengths, these values are still much lower than those seen in experiments, with typical run lengths of 800 nm and 1.5 mm measured for murine and bovine dynein in vitro [69]. This suggests that although some processivity can be achieved with independent head domains, coordination is  important to obtain the higher run lengths that are seen in experiments. The fact that the mean run lengths are lower than those seen in experiments could be due to the regulatory functions of dynactin and other cargo adaptor proteins such as BICD2 present in vivo, which activate long-distance movement of the motor [33,34,42]. Further modelling in this direction might help to confirm or refute such hypotheses. Variations in g ATP have little effect on the percentage of not passing and alternating steps; however, they do have an effect on run length, with an increase in g ATP . 10 MDa ns 21 leading to a fall in the mean number of steps and lower run length for all variables (tables 5 and 6).
Remark 5.4. We have extended this model framework to include random backward stepping and a variable step size, details are discussed below in §5.4. If independent stepping is assumed, then this leads to a reduction in run length to 158 nm for the cargo. This is likely to be due to the presence of backward steps shortening the run length. However, larger step sizes may also lead to an increase in detachment time for a single head within the model, increasing the likelihood that the other head will also detach. MTBDs D and E   In this section, we extend our modelling framework to take into account the backward stepping, variable step sizes and large-scale dwelling of dynein. Previously, we used a fixed time interval T Final /N for the stepping of a single MTBD; however, the active stepping of the MTBD should end when the MTBD binds to the microtubule. Consider the interval where p kþ2 is the next binding site for the unbound MTBD j and t max is the maximum potential length of the stepping interval. Hence, the total time spent stepping is given by . In order to model dwelling over large timescales, we take a multiscale approach by using one timescale for stepping and one for dwelling; variable step sizes and backward stepping are also included in the model (see the electronic supplementary material for details).
The model is solved numerically in MATLAB for N ¼ 100 steps using ode45 for the stepping model and the stiff solver ode15s [68] for the dwelling model. The non-dimensional systems are solved and then converted back to dimensional results so that they can be presented together on the same timescale. Table 5. Mean percentage of not passing and alternating steps, and mean number of steps in a run given a range of values for the parameters K ATP (MDa ns 22 ) and g ATP (MDa ns 21 ). The data represent the results of 100 simulations with a mean dwell time of m ¼ 2 Â 10 9 ns. If x% of steps are not passing, then (100 2 x)% of steps are passing. Similarly, if x% of steps are alternating, then (100 2 x)% of steps are not alternating; except for the case labelled* where the number of steps in a run was always less than or equal to one and hence neither alternating or non-alternating steps were present.  The maximum length of the stepping time interval is taken to be t max ¼ 10 6 . The primary values are taken as follows: the probability that MTBD E steps given that MTBD D (E) stepped previously is taken to be P D ¼ 84% (P E ¼ 16%), the maximum separation distance is taken to be 48 nm and the probability of backward stepping (P Back ) is set to be 20%. The mean dwell time is taken to be m ¼ 2 Â 10 9 ns as experimental results have shown the average dwell time for dynein to be 2 s [48]. For each step, we take n from the Poisson distribution about 2 to give the step size nL ATP for the forward step sizes and n from the Poisson distribution about 1 for the backward steps to give the step size 2nL ATP (see remark 5.5). We assume that zero steps are possible, but they are not counted toward alternating or non-alternating steps. See table 2 for the range of values or distributions used for the stochastic parameters.
Remark 5.5. Note that the distribution used to obtain the step size nL ATP could be obtained through analysis of the experimental data to give a more accurate representation of the step sizes of a particular dynein species. However, it could also be used to analyse the effect of different distributions on stepping behaviour and run lengths which is left for future studies.
The results show similar profiles for the tail, AAAþ rings and MTBDs, with a clear presence of backward steps, variable step sizes and increased dwell times between steps; however, we see a significant difference for the velocity of the cargo (figure 11). Our computational results show that the frequency of alternating steps does not differ and is an emergent process of the modelling. On the other hand, our results show that the not-passing steps differed by approximately 1-4% for optimal parameters, which is not such a big variation. However, this variation becomes significant for small values of the maximum separation distance (see the electronic supplementary material for details). The cargo now dwells between steps, with an oscillatory velocity profile that returns to zero between steps which is similar to the in vivo experimental results shown by Garrett et al. [10]. Using the primary values for the parameters gives a maximum velocity of the cargo of 15 Â 10 5 nm s 21 , and a velocity of up to 2 Â 10 8 nm s 21 , for the tail domain. This is much higher than velocities measured experimentally with dynein typically moving at speeds of 600 nm s 21 , at saturating ATP levels and at room temperature with in vivo velocities reaching up to 3 mm s 21 in mammalian neurons, although yeast dynein moves at slower speeds of around 50-80 nm s 21 [69]. A full parameter analysis of all unknown model parameters needs to be conducted in order to establish the parameter set which gives quantitatively accurate values for the velocity for each species and context.
Although the overall direction of travel for the AAAþ rings and MTBDs are closely related, we do see differences in their behaviour, with the two AAAþ rings being further apart from each other than the two MTBDs and crossing paths at different time points (figure 11e,f ). This would suggest that Table 6. Mean run lengths for the cargo and tail domain given a range of values for the parameters K ATP (MDa ns 22 ) and g ATP (MDa ns 21 ). The data represent the results of 100 simulations with a mean dwell time of m ¼ 2 Â 10 9 ns. If x% of steps are not passing, then (100 2 x)% of steps are passing. Similarly, if x% of steps are alternating, then (100 2 x)% of steps are not alternating. We also achieve a range of backward steps in the model. For a fixed step size and minimal dwelling of 2 ns (see remark 5.6), taking the probability of backward stepping to be 20% and the maximum separation distance to be 48 nm results in 24.91% backward steps, reflecting the experimental observations of Qiu et al. (23%) [37]. While in order to match the observations of Reck-Peterson et al. (13%), we can take the probability of backward stepping to be 10% and the maximum separation distance to be 56 nm to obtain 13.19% backward steps [48]. A probability of backward stepping of 0% does not mean that there will be no backward stepping in the model as we use this parameter to represent random backward stepping, which we differentiate from the corrective backward steps taken when the MTBDs are too far apart. We see from the results that this gives a very low presence of backward stepping, much lower than in experimental results. Hence, this suggests that the MTBDs might randomly step backwards or that tension within the complex causes restorative backward steps through some mechanism not explicitly modelled. However, such processes are beyond the scope of this study. Experimental studies have shown that dynactin plays an important role in the directionality of dynein, and hence we may need to explore these effects in greater detail [2,6,[28][29][30]32,34,65]. The effects of stepping along the microtubule in two dimensions may play a role in backward stepping if the motor domain rotates due to the off-axis components of the steps. In our current model, we are setting external forces to be zero, but these forces may play a role in the directionality of the head domain for in vivo studies.
Remark 5.6. The effects of backward stepping and variable step sizes were explored on an initial minimal dwelling model for a single timescale, using the non-dimensionalization given in §4.2, taking m ¼ 2 ns.
We explored variations in the maximum separation distance on the stepping patterns (table 7). The reduction in maximum separation distance increases the likelihood of backward stepping, this is to be expected as backward stepping is directly related to the separation distance in the model, with an unbound head stepping backwards if it is too far in front of the other MTBD. We also see that reducing the maximum separation distance increases the likelihood of passing steps, which makes sense as closer MTBDs are more likely to cross over one another during stepping.
Increasing the stepping probability of MTBD E after MTBD D has stepped decreases the percentage of not passing steps and also decreases the percentage of backward steps (table 8). This is likely to occur as the increased coordination would create a more efficient stepping pattern reducing the prevalence of wasteful backward steps by keeping the motor domains closer together and hence passing steps would also be more likely to occur. Owing to the presence of zero sized steps, in order to achieve similar results to experiments we take the probability that MTBD E steps set at 84% if the previous step was taken by MTBD D, and 16% otherwise. This results in 82.72% non-passing steps, 74.68% alternating steps and 20.91% backward steps ( figure 12). Approximately 10% of steps by the MTBDs were of a zero step size, the majority of steps were of 8 -16 nm and histograms of both forward and backward step distributions (not including zero steps) are given in figure 13. It is likely that dynein does experience step sizes of 'zero' length, i.e. detaches but rebinds to the same point on the microtubule. However, this stepping behaviour is not picked up (and therefore not accounted for) by the step-finding algorithms used in experimental data analysis. This suggests that the coordination between head domains could actually be higher in reality than recorded in experiments.

Discussion
In this study, we have derived a general integrative mechanistic model that describes the transport mechanism of cytoplasmic dynein. Our results give a mixed stepping pattern with a predominantly not passing stepping profile, which is an emergent process of the model, matching experimental Table 7. Mean percentage of not passing, alternating and backward steps given a range of values for the maximum separation distance d (nm). The data represent the results of 100 simulations with the probability that MTBD E steps set at 74% if MTBD D stepped previously and 26% otherwise. The probability of random backward stepping is set to be 10% and the mean dwell time is taken to be 2 ns. If x% of steps are not passing, then (100 2 x)% of steps are passing. Similarly, if x% of steps are alternating, then (100 2 x)% of steps are not alternating. observations. We have shown that there is likely to be some form of interhead coordination between the timing of the ATP cycles of the two AAAþ rings in order to account for the alternating patterns seen in experiments.
We have been able to model uncoordinated motion and have shown that dynein can still achieve some level of processivity through this mechanism. For example, the model achieves run lengths close to those seen for murine dynein in the absence of dynactin, when using a fixed forward stepping pattern. This may suggest that either dynactin has some influence on the coordination of the motor domains or that we need to account for the effect of dynactin in our model in some other way.
Loa dynein in mice has been shown to exhibit shorter run lengths than wild-type complexes [10][11][12]. Ori-McKenney et al. [11] measured run lengths of 259 nm for Loaþ/2 mutants and 175 nm for Loa2/2 mutants; we are able to achieve similar run lengths through an appropriate choice of parameters. Ori-McKenney et al. [11] suggest that the Loa mutation may cause altered coordination in the motor domain of dynein. Our results suggest that it may be possible that this mutation disrupts the coordination within the complex, potentially leading to more frequent detachment of the motor from the microtubule and shorter run lengths. Deng et al. [14] have also shown that the Loa mutation causes dynein to have a lower affinity to dynactin, and so it may be through this disruption that the mutation affects the transport mechanisms of dynein.
It would be interesting to investigate the effect of the differences in these dwell times rather than assuming that the mean dwell time for each head domain is equivalent. In particular, allowing the lagging and leading head to have different dwell times may encourage a more coordinated stepping pattern, and experiments have shown that the lagging and leading heads have different stepping characteristics [36,37].
Currently, the motor domain can diverge as we assume that once the MTBD is bound, it is bound until a conformational change through ATP hydrolysis cycle occurs, it would therefore be interesting to introduce the effect of forces on detachment in this model. It has been shown by Gennerich et al. if MTBD D had stepped previously and 16% otherwise. The maximum separation distance is set to be 48 nm, the mean dwell time is 2 ns and the probability of random backward stepping is taken to be 20%. The maximum separation distance is set to be 48 nm, the mean dwell time is 2 ns and the probability of random backward stepping is taken to be 20%.
rsos.royalsocietypublishing.org R. Soc. open sci. 5: 171568 [47] that dynein can walk through applied force alone, and so these forces are important to the model. Another alternative is to introduce explicit modelling of the linker that has been shown to gate the ATPdependent release of dynein from microtubules [4]. In particular, the linker has been shown to play a critical role in gating dynein stepping behaviour at high interhead separation distances. By introducing a linker, a tension-dependent force that acts to retract the leading head or to pull the lagging head will counterbalance the larger separation distances between heads during dynein stepping behaviour. We have also been able to incorporate backward stepping, a variable step size and dwelling over large timescales into the model. The results give trajectories for the complex and cargo that qualitatively match experimental observations. Although we have compared our results qualitatively and quantitatively with results published in the literature, it would be beneficial to carry out detailed comparisons for a specific dynein whereby space-time series distributions data are provided. An ideal candidate is to employ a Bayesian parameter identification approach that allows us to compute optimal parameter distributions resulting from fitting the solution of the mathematical model (with all parameters assumed unknown) to experimental data (known) in an optimal sense [70]. The result of this approach is the rich statistical data that provide various statistical measures such as mean, variance and 95% credible regions. Furthermore, velocities can also be computed as distributions, which is more suitable for analysis and comparison to experiments. This approach forms part of our current studies, the only requirement is finding appropriate experimental data generated in terms of space-time series to allow us to optimize parameter identification such that the model solution best fits the data.
We have also shown that backward stepping that is directly related to the separation within the complex cannot account for the high percentages of backward stepping seen experimentally, and hence there must be something else external to this simple model causing these characteristics. We suggest that the impact of dynactin on the transport mechanisms and the three-dimensional nature of dynein need to be explored further with regard to their impact on backward stepping.
By prescribing the levels of coordination within the model, we can match experimental observations of the alternating stepping pattern, but when considering the possibility of 'zero' step sizes, this coordination must be higher than that seen in experimental observations. The model predicts the preference of dynein to the not passing stepping pattern when the motor is allowed to separate (which is realistic due to the large step sizes seen in experiments), and this matches experimental observations. The model also predicts that species of dynein which prefer a tighter conformation may be more likely to experience backward steps and have a higher prevalence of passing steps. Stronger coordination between the two motor domains could also reduce backward stepping, which leads to more efficient stepping as backward steps may be wasteful.
Other future works involve studies to investigate the effects of variable dwell times and straindependent stepping to establish a complete model which incorporates all aspects of transport mechanisms for cytoplasmic dynein. Apart from extending the model to take into account dynactin and tension linker domain, it will also be interesting to model multiple dyneins and how they aid or hinder the stepping behaviour. Recent studies by Urnavicius et al. [28] and Grotjahn et al. [29] reveal that dynactin has the capacity to recruit a team of dyneins for processive motility. However, to the best of our knowledge, no mathematical model has been formulated that could describe such experimental observations. The approach presented here sets foundations for such a study.