Moving in the mesoscale: Understanding the mechanics of cytoskeletal molecular motors by combining mesoscale simulations with imaging

Rapid advances in experimental biophysical techniques are generating a wealth of information about the mechanical operation of the cellular cytoskeleton and its motors. However, each of these tools typically provides only a limited piece of a highly complex puzzle. There is a need to develop new computational tools that can integrate these data together into a central model. Here we discuss the experimental advances alongside the computational tools, and propose how these could be developed to successfully combine the emerging structural and dynamic experimental data on cytoskeletal motors. We consider examples of both single motors and arrays of motors within a biological cell.

functioning cell. While much of our information about the cellular scale organization of the cytoskeleton has come from light microscopy, molecular details of motor proteins have come predominantly from X-ray crystallography or nuclear magnetic resonance, and more recently cryo-electron microscopy (cryoEM). Consequently, there has been a gap in experimental resolution between length scales of 10 and 500 nm and 1 μs-1 s, which is the regime that we define here as the biological mesoscale ( Figure 2a).
The cytoskeleton is made up of microtubule (MT) and actin filaments (F-actin), which provide the tracks for cytoskeletal motors to work along. Kinesin and dynein operate along MTs while myosin uses F-actin. 1 Motors function through diverse mechanisms, either as part of highly organized super-macromolecular arrays, 2 or on an individual/ more random group basis. 3,4 For example, axonemal dyneins assemble into highly ordered arrays that connect adjacent MTs in the axoneme of motile flagellar and sensory cilia. Their synchronized motion results in bending of the axoneme to drive ciliary and flagellar beating ( Figure 1). [5][6][7] Myo5a however works independently or as part of a group to transport molecular cargoes to the edges of the cell, 4 such as melanosomes in skin cells 8 (Figure 1). It is a processive motor that "walks" in a hand-over-hand manner along its F-actin track. 9

| EXPERIMENTAL METHODS FOR THE STUDY OF CYTOSKELETAL MOTORS
There are a growing number of biophysical techniques, with ever increasing spatiotemporal resolution, which provide new insight into the mechanics of cytoskeletal motors (and large biomolecular complexes in general), the cellular environments in which they need to operate, and their full range of biological functions ( Figure 2). Microscopy and force spectroscopy have played an important role in visualizing motor proteins and generating hypotheses on the mechanics and coordination of their motion. Single molecule techniques enable targeted analysis of individual molecules, this is important for gaining an understanding of the biomolecular mechanics through visualization and quantification of biophysical properties, such as binding forces. 10

| Electron microscopy and tomography
Electron microscopy (EM) enables visualization of proteins trapped in different dynamic states. It provides structural information but cannot probe rates of change during biomolecular processes (although some rates can be inferred from the population of conformations occupied). Negative stain EM (nsEM) has historically provided a wide range of information on conformational states of motor proteins, 20 insight into how motors interact with their tracks, 21 mechanical properties, 22,23 and even binding kinetics. 24 However, the structural resolution of nsEM is limited to $10 Å, 25 while that F I G U R E 1 Schematic representation of a super-macromolecular array of dyneins within a primary cilium (left), and myosin 5a (Myo5a) working independently to move a cargo along the actin network (right) of cryoEM (Figure 2c,d) is now beginning to rival X-ray crystallography (currently 1.22 Å best resolution 14 ) as a result of better detectors and advanced processing techniques. 26 CryoEM has now been used successfully to resolve the structure of a full-length motor protein for the first time, 27 motors bound to their tracks, [28][29][30] and the organization of actin and myosin in muscle myofibrils, 31 none of which would have been possible by X-ray crystallography. CryoEM structures may provide a closer insight to the native state of the protein as the molecules are preserved in an aqueous environment, 32 and fleximers are able to move freely up to the point of freezing the specimen, unrestricted by crystal packing and artefactual interactions. [33][34][35] Time-resolved cryoEM can capture dynamic states at different time points during a reaction to determine the changes in structure during this sequence of events. 36,37 Cryo-electron tomography (cryoET) allows in situ visualization of proteins in their native environment, through sectioning and imaging regions of interest directly from cells, revealing the complexity of the environment through which motors can move albeit at lower resolution than cryoEM. 38,39 2.2 | Light microscopy Light microscopy, in particular "super-resolution" techniques, provides information about the position, abundance, and organization of specific labeled proteins, but only low-resolution information about biomolecular structures (Figure 2h). Many in vitro experiments reconstitute the motors and the tracks, but do not take place in cellular conditions or include F I G U R E 2 An indication of scale. (a) Scale describing the temporal and spatial resolution of experimental techniques and the "biological mesoscale" (10-500 nm, 1 μs-1 s). Highest temporal resolution of super-resolution microscopy = 50 μs (interferometric light scattering microscopy 11 ), and force spectroscopy = 67 ms (high speed atomic force microscopy 12 ). Highest spatial resolution of superresolution microscopy = 1 nm (single molecular light microscopy 13 ), and electron microscopy (EM) = 1.2 Å (cryoEM 14 ). Arrowless heads indicate the resolution range goes beyond the scale displayed. (b) "Bead and spring" elastic network model (ENM) of Myo5a motor domain (PDB: 1OE9 15 ) using ProDy. 16,17 Arrows show anisotropic network model modes of motion. Dotted lines indicate division of motions as domains and larger beads show "hinge" regions motions predicted from Gaussian network model. Yellow = mode 1, green = mode 2, pink = mode 3. Networks (MEDYAN) model of actin filaments bundling after 2000 s due to contractile forces from motors and cross-linkers. Simulation box volume is 1 Â 1 Â 1 μm 3 . Blue = actin filament, orange = myosin 2 mini-filaments, green = α-actinin cross-linker. Data displayed from. 19 (h) Confocal image of F-actin in a cell cargo. Light microscopy of intact cells can, however, take place in a native environment and can determine how the motors generate a specific biological outcome, such as cell trafficking or construction (e.g., of the axoneme). 40,41 Fluorescence microscopy has also been used to determine the impact of track modifications on motor function. 40 There is now a variety of techniques for imaging the positions, relationships, and dynamics of multiple targets within cells. "Super-resolution" information can now be obtained well beyond the classical diffraction limit of a microscope ($250 nm), down to a length scale of $1 nm 13,42 (Figure 2a). However, there is often a trade-off between spatial and temporal resolution that has driven the development of systems that capture fine detail at higher speed. [43][44][45][46][47][48] Where molecules are found in ordered arrays, single molecule fluorescence localization microscopy can be analyzed to extract their arrangement. 49,50

| Single molecule techniques
Although single molecule techniques often do not take place in a native environment, they have been instrumental in visualizing molecular dynamics (MD) and elucidating physical properties of motor proteins. In optical microscopy, the use of total internal reflection fluorescence microscopy (TIRFM) in particular allows dynamic imaging of labeled single molecules with a precision of $10 nm at $100 frames per second. 51 Interferometric scattering microscopy can provide even higher spatiotemporal resolution (5 nm and 50 μs) 52 ( Figure 2a).
Other methods probe the interactions and environment of molecules. Single molecule force spectroscopy methods include the use of atomic force microscopy (AFM) and optical and magnetic traps, and probe the forces between binding partners. [53][54][55] Förster resonance energy transfer can also be used to measure short distances (<$10 nm) between molecules, or two parts of the same molecule, enabling force sensing as this distance changes. 56 Connecting structural (e.g., cryoEM) and physical information (e.g., optical trap) to dynamics (e.g., TIRFM) and understanding in the full cellular context (e.g., super-resolution light microscopy) requires an understanding of biomolecular mechanics.

| MODELING METHODS FOR CYTOSKELETAL MOTORS
Two distinct computational approaches have been used to model cytoskeletal molecular motors; spatially resolved models, where proteins are represented explicitly as physical objects, and kinetic or reaction-diffusion based models, where differential equations that describe chemical rates are used to evolve the system as a function of time. While here we concentrate specifically on spatially resolved simulations, much insight has also been obtained from simulations that represent processes (rather than objects) using sets of kinetic equations 57,58 and in some cases, these have been coupled to spatially resolved models (e.g., MEDYAN model of the cytoskeleton, Figure 2g, described below).
The importance of multiscale simulations of the cytoskeleton and the diversity of computational tools needed to span the atomistic up to cellular length scales is emphasized in a number of recent reviews. Focusing on dynein, Dutta et al. provide an overview of the lessons learned from computer models, spanning quantum mechanical calculations probing the catalytic mechanism, up to length scales where a dynein dimer is modeled on an MT track. 59 Recent summaries of the simulation tools available for mesoscale simulations, covering cytoskeletal simulations, protein aggregation, and nuclear dynamics are available, 60,61 the computational tools for modeling whole cell cytoskeletal mechanics have also been reviewed. 62 3.1 | Atomistic and coarse-grained particle models In biomolecular simulations, particles can be used to represent single atoms (e.g., atomistic MD simulations), groups of atoms such as amino acids, or larger biomolecular assembles such as tubulin dimers 63 and segments of actin filaments. 19 MD simulations use a bespoke force-field to capture the complex chemical interactions between every atom within the protein, which are parameterized from the "bottom-up" using a combination of quantum chemical calculations and experimental (e.g., vibrational spectroscopy) data. In principle, atomistic MD simulations contain sufficient information to capture the sequence-dependent mechanics of motor proteins, such as conformational changes associated with track binding, ATP binding and hydrolysis, and the explicit effect of the macromolecular environment. However, the computational expense of the calculations limits the size of molecule and timescales that can be explored (typically 10 5 atoms and $1 μs), so that currently observing such changes in an atomistic simulation is impossible, except with advanced sampling tools. For example, Fischer et al. use conjugate peak refinement to calculate a minimum energy path to enable simulation of the recovery stroke in myosin (1-10 ms timescale). 64 Additional insight into features of the power and recovery stroke in the motor domain of myosins, [65][66][67][68] and the atomic details of track-motor interactions between actin and myosin, 69 has been obtained from MD simulations. Moreover, they provide a means of fitting models to electron density maps in a way that enforces chemical laws, which is particularly useful for low resolution maps often seen in flexible proteins. [70][71][72] A combination of atomistic and coarse-grained MD simulations (with Cα atoms represented only) has been used to understand the effect of posttranslational modifications on the diffusion of proteins on the surface of MTs. 73 The atomistic simulations were used to parameterize coarse-grained models of the wildtype and polyglutamate-modified tubulin tails. The use of coarse-grained representation allowed simulations of a 4 Â 6 MT lattice, which would have been prohibitively expensive at the atomic level. The simulations revealed the crucial role of tubulin tails in mediating surface interactions; specifically those that effect electrostatic interactions.
Coarse-grained elastic network models (ENMs) typically have Cα atoms represented as beads and interatomic interactions within the protein represented as simple harmonic springs to facilitate the solution of the equations of motion ( Figure 2b). They are used to model slow-mode dynamics (μs to ms timescale 74 ) and are computationally cheap compared to atomistic simulations. The spring constant and connectivity can be chosen to represent internal protein structure, and so can capture key aspects of sequence dependent mechanics, such as the positions of flexible hinges within the biomolecule. 75 This approach has been used previously for Myo5a to provide a structural explanation of the impact of strain on motor conformation, 76 the mechanism by which ATP binding leads to opening of the motor cleft, 77 and the kinetic coordination of the dimer. 78 For cytoplasmic dynein, a multiscale scheme combining atomistic, coarse-grained, and phenomenological models has been devised to show how tension between the two dynein heads modulates track binding efficiency. In this model, the tension causes the trailing head to detach more easily from the MT, while the leading head remains bound, resulting in a forward stepping motion. 79 ENMs can also be used to fit atomistic models into cryoEM maps, 70,72,80 and can now be generated from cryoEM maps to predict even coarser conformational dynamics. 81 The dynamic behavior of MT tips has been captured by representing each tubulin monomer as a particle, and then arranging tubulin dimers into a cylindrical lattice. 63 While the bonds between monomers in a dimer were unbreakable (e.g., harmonic springs), interaction potentials between dimers allowed for disassociation while new dimers were kinetically associated to the tip, which is essential for modeling MT growth and catastrophe. By defining different potentials between GTP and GDP tubulin within the MT lattice, it was possible to study the stabilizing properties of the GTP cap. This study demonstrates that active processes that consume chemical energy can be modeled at the coarse-grained level, without the need to explicitly represent protein conformational changes.
Ultra-coarse-grained particle models have also been used to provide insight into dynamics on the larger spatiotemporal scale necessary to model cytoskeletal reorganization in the context of the whole cell or specific cellular architectures, such as cellular protrusions. AFINES simulations explore timescales of seconds and tens of μm length scales, 82 F-actin is represented as 2D chains of beads connected by springs and neglects steric interactions, limiting its application to systems with low volume fractions of F-actin. By adjusting filament lengths and concentrations of contractile elements, simulations of actomyosin networks revealed three structural phases that were bundled, polarity-sorted, and contracted. 83

| Continuum mechanics models of cytoskeletal motors
Continuum mechanics simulations of cytoskeletal motors represent the complex 3D geometries of proteins or larger subcellular structures as a three-dimensional continuum, typically subdivided into simple domains or elements using a mesh. Long thin objects, which are prevalent in the cytoskeleton (e.g., MTs, actin), can represented as sets of onedimensional rods. These meshes/rods are then used to compute the mechanics of the object, by solving the equations of motion for the system for each node within the mesh or rod.
In continuum mechanics internal forces are determined by the constitutive model, which relates the forces in the system to its deformation. Commonly, an elastic constitutive model is used which, similar to an ENM, describes a system that will return to an equilibrium structure if external forces are removed. It is straightforward to add simple viscous dissipation to this description. Yet, a huge array of constitutive models are available, including effects such as viscoelasticity for materials with time dependent memory of deformation, or plasticity for materials which flow only above a critical force. Analogous descriptions apply to lower dimensional models such as 1D rods, which commonly include elastic forces for stretching, bending, or twisting and viscous resistance to motion 84 but rarely go beyond that.
At the lower length scales accessible by continuum mechanics, where thermal fluctuations are significant, the fluctuating finite element analysis (FFEA) method aims to represent both the mechanics of motors, and their changing shape due to thermal fluctuations. 85,86 The equilibrium molecular shape is obtained by cryoEM or cryoET, converted into a 3D finite element (FE) mesh, and/or assembly of 1D rods for elongated biomolecules (Figure 2e). 84 A simple viscoelastic constitutive model is used to represent the material properties, which can be tuned to capture specific sequence selective properties, such as hinges within a stiffer structure. However, unlike atomistic simulations, or ENM models, these properties need to be input into the model as parameters, obtained either from more detailed atomistic simulations, 87 or experimental data. 88 Interactive forces, such as protein-protein interactions, are represented by short range (e.g., Lennard-Jones) potentials and must also be input as parameters. Conformational changes are modeled throughout simulations in FFEA by interpolating between meshes representing the distinct equilibrium structures of a chemo-mechanical cycle. 89 FFEA has been used to study a variety of molecular motors, including both cytoplasmic 87 and axonemal dynein, 89 where simulations mimicking several iterations of the chemo-mechanical cycle could be performed.
At an equivalent length scale to the particle-based AFINES software, the 3D MEDYAN 19 ( Figure 2g) and 1/2/3D Cytosim 90,91 software packages represent cytoskeletal filaments as flexible chains of interconnected rods with finite and infinitesimal widths, respectively, that can sterically interact. The dynamics of filaments and other diffusing objects such as molecular motors are governed in Cytosim and AFINES by overdamped Langevin dynamics, and in MEDYAN by a stochastic reaction-diffusion scheme coupled to mechanical energy minimization. Across all three models, the mechanical bonds formed by motors and cross-linkers between adjacent filaments are represented implicitly as harmonic potentials with cross-linkers being stationary and motors dynamic once they are bound. Virtual AFM simulations with MEDYAN revealed actin bundling and alignment in the direction of applied force was mediated by both mechanical and chemical cytoskeletal activity at fast and slow timescales, respectively. 92 A recent application of Cytosim investigated how the coarse graining of different ensembles of non-muscle Myosin 2 affected their contractile dynamics. 93 Whole cell models have been used to investigate the role of actin and myosin in cell motility. A whole cell was meshed in two dimensions with actin and myosin represented kinetically to model cellular migration, 94 and analytical flow fields solved inside a deformable sphere were able to predict contractile flows observed in Zebrafish oocytes. 95 A review of the computational techniques used to understand how cellular mechanics arises from the underlying cytoskeleton is available. 62

| COMBINING EXPERIMENT AND SIMULATION: WHAT PHYSICS IS NEEDED?
The experimental and computational tools that have been applied to cytoskeletal motors have provided numerous pieces of a very complex puzzle (e.g., structures/molecular abundance/forces/kinetic rates) across a variety of spatial and temporal resolutions (e.g., Å to mm and μs to s), as shown in Figure 2 and discussed in Sections 2 and 3. A comprehensive understanding can therefore only be obtained by combining complementary biophysical data together. This requires a common model for the mechanics of motors that captures the physical processes underlying the data obtained, which can only be achieved using multiscale modeling techniques. An integrative model is most useful when it focuses on the length and/or timescales where distinct experimental techniques overlap. For example, the length and timescales at which super-resolution microscopy and cryoEM/ET can be connected are in the range of 10 nm (for a single motor) up to μm for processions of motors along tracks.
At any given resolution, the computational expense depends critically upon the physics needed to capture the biomolecular mechanics of interest. ENMs are extremely efficient, because the particular "harmonic balls and springs" form of the equations used in these models makes the dynamics fast to compute, even at atomic resolution. However, this is not the case for general potentials that are not simple harmonic springs. Moreover, if electrostatics or hydrodynamics are represented explicitly, then the computational expense is massively increased, because these are long range interactions that necessitate the consideration of large volumes of the simulation cell. Simulations requiring long range interactions are also more difficult to parallelize. These considerations emphasize the importance of establishing the appropriate physics of the system.
We now discuss the physics needed for simulations of cytoskeletal motors at various length scales. We also consider how these representations could be combined into a holistic multiscale representation from the atomistic level through to the cell.

| Nature of the irreducible model elements
The construction of a hierarchical multiscale scheme for representing cytoskeletal motors requires the nature of the "irreducible elements" at each length scale to be defined, which in turn dictates the physics needed. For example, models of catastrophe and regrowth of individual MTs 63 represent α and β-tubulin protein subunits in the MTs explicitly, but at a more coarse-grained level the dynamic nature of F-actin and MTs are represented by extensible 19 or inextensible 90 elastic rods. The latter provides less detailed information, but accesses biological processes at larger length and timescales, such as cytoskeletal reorganization.
The level of representation determines the definition of "solute" and "solvent" within the simulations. Options range from a fully atomistic representation of all of the waters, counter-ions, metabolites (e.g., urea), and bioactive molecules such as drugs and hormones, through to an entirely continuum representation, in which the whole external environment is treated implicitly. The physics needed to describe the environment depends on the length scale represented. For example, while it is crucial to accurately represent electrostatic interactions in atomistic simulations, at the mesoscale these may be effectively screened by mobile counterions, allowing computationally cheaper, short-ranged potentials to be used. At larger length scales, the environment may be the crowded cytoplasm, which will contain a complex mixture of proteins and metabolites. 96 To construct a simulation of a representative section of the cytoplasm in a living cell, where numerous cytoskeletal molecular motors are continuously active outside of the region of interest, the cumulative effect of these external active processes may need to be included. Force-spectrum-microscopy performed in living cells has shown that force fluctuations from molecular motors can enhance the motion of cellular components, and has revealed that at the cellular level the cytoplasm behaves as a weak elastic gel, which can the drive fluctuations in objects embedded in this network with characteristic timescales that depend on the operational timescales of the motors themselves. 97 For whole cell models using continuum FE models the cytoplasm and cytoskeleton combined will be characterized simply by stiffness, viscosity, the choice of constitutive model, and the resolution of the mesh. 62 One major question for coarse-grained simulations is how to represent intrinsically disordered proteins. These regions have important biological functions over multiple length scales that are poorly understood. From the viewpoint of protein sequence and chemical structure, disordered regions provide a diffusive tether connecting protein subunits. From a thermodynamic view focused on subcellular length scales, disordered tails display fascinating emergent behavior that has been hypothesized to drive liquid-liquid phase separation, and to establish assemblies of "nano-factories" of biomolecular machines that are localized in membrane-less organelles through entropic effects. 98 Short, disordered tails decorate the exterior of MTs, whereas in other biomolecular machinery, including the kinetochore, these tails can be extremely long and so occupy significant volume within a crowded cellular environment. 99,100 Exploring the role of disordered tails in mediating mesoscale behavior in ultra-coarse-grained biomolecular simulations is likely to reveal exciting new physics of direct relevance to biological function.
As well as reducing the computational expense, models that use a less detailed representation also produce smaller datasets, which are far easier to store, visualize, and analyze. Atomistic MD trajectories in explicit water can be terabytes in size, and require considerable extra computational effort to process before meaningful insight can be obtained. For biological structures at the mesoscale, such as a collection of MTs and their cytoskeletal motors, these data sizes are currently prohibitively large for atomistic simulations.

| Molecular shape and ultrastructure
Spatially resolved models must provide a representation of the objects they are describing, which requires knowledge of the molecular shapes and organization at the level of interest. In the cell, molecular motors act in coordinated super-macromolecular arrays (e.g., axonemal dynein and muscle myosins), or as collectives of molecules carrying large cargos (e.g., kinesin, cytoplasmic dynein, Myo5a). They may contain variable numbers of motors in any given motor complex, and cofactors such as dynactin may also be required. 101 CryoEM/ET can provide information about the super-macromolecular arrangements of proteins into large protein complexes, and how groups of motors coordinate transport. However, it is not always straightforward to identify the components in the density maps, and abstract them into a physical model, particularly if only low-resolution information is available. Complementary information about the positions, abundances, and coordinated motion of particular proteins can be obtained from super-resolution light microscopy studies. 102 Computational tools for identifying species in cryoET maps will massively assist in model building. For macromolecular complexes where atomistic information is available, flexible fitting of atomic coordinates into the cryoET densities opens up the possibility of mesoscale and ultrastructural models that are informed by atomistic information, and potentially coarse-grained models at larger length scales. 81 Automated procedures and workflows for constructing models will greatly facilitate this approach (Section 5.1).

| Structure dependent biomolecular mechanics
Individual biomolecules are soft, nanoscale objects that change shape significantly due to thermal fluctuations. The nature and magnitude of these conformational changes are dictated by the underlying atomistic structure of the protein. Even in the absence of energy consumption, thermal motion can involve global conformational changes across the whole molecule. Atomistic structure impinges upon the mesoscale whenever it influences mechanics through changes in secondary structure, for example if a helix region transitions into a disordered loop, or if secondary structural elements change their relative positions within the protein. Communication of ATP driven structural changes is essential to motor function, for example in myosin motors ATP hydrolysis is associated with kinking of a helix (the "relay helix"), partial closure of the actin binding cleft, and weak actin binding. 103 While, P i release is associated with straightening of the relay helix, twisting of a β-sheet, further closure of the actin binding cleft, and strong actin binding affinity. 30 In turn, the stiffness, organization and interactions between mesoscale cytoskeletal components determine cellular mechanics, for example, the viscoelastic responses of different cellular regions.
The mechanical properties needed to model actin and/or MT assemblies have been measured experimentally, because these filaments are large and stable enough to be visualized with AFM, negative-stain, and cryoEM, and have also been studied using complementary biophysical tools, including single molecule techniques. [104][105][106][107][108][109][110] However, for smaller cytoskeletal structures, and for molecular motor components, obtaining suitable parameters for modeling biomolecular flexibility is more challenging. Although, with novel cryoEM image processing techniques 111 that determine molecular flexibility, this may become easier.
Biophysical investigations into the walking mechanism of motors have shown that these parameters are key to function. For example, in Myo5a, lever stiffness is a key part of its mechanism. 112 Bending of the lever in the leading head when both heads are bound to F-actin can be seen in nsEM 21,24 and cryoEM ( Figure 2d) providing a visual indication of the intramolecular strain that enables processive walking. Insight into the specific physics of this strain has been gained using single molecule techniques that measure stall and binding forces directly, such as optical trap and AFM. [113][114][115][116][117] Strain also biases ADP dissociation kinetics in each head, accelerating release in the trail and slowing release in the lead head. 24,118,119 Additionally, intramolecular strain has been shown to be responsible for the lever swinging motion required for walking, 52,114 and not chemical energy liberated from ATP hydrolysis, via AFM. 120 Similarly, multi-resolution simulations of cytoplasmic dynein have shown that the directional response of the dynein stalk to applied tension is key to the mechanism of coordination of the two dynein heads during stepping. 79 Based on these experimental observations, a spatially resolved model of molecular motor activity, such as the walking of cytoplasmic dynein or Myo5a, needs to capture both the tension stored by the different molecular configurations within the chemo-mechanical cycle, and the flexibilities (including directionality) of distinct molecular hinges. When atomistic information is available, hinge regions can be identified by locating unstructured regions within the secondary structure of the protein, for example, proline residues break up the coiled-coil registry in the stalk of dynein. 121 The simple shape of the protein, for example, whether it is globular, or long and thin, also plays a significant role, which is why ENMs provide such a successful coarse-grained representation of protein dynamics.
For particle based coarse-grained models, including ENMs, this level of detail can be captured by carefully choosing the connectivity, and potentials between particles. In continuum FE, FFEA, and rod models, sequence dependent mechanics can be captured by assigning different material parameters local to a region to define the inhomogeneous (and direction dependent) mechanics of the protein. However, this requires a separate parameterization procedure, either by running atomistic MD simulations on tractable fragments, or by comparison with experimental data, such as nsEM.
Given the enormous number of proteins in eukaryotic cells, the parameterization of mesoscale coarse-grained models needs to be automated. In ENMs, the dynamics depends on how many pseudo-particles are used to represent the protein, and how the harmonic network is connected. For proteins that have a known atomistic structure, this is in principle assignable algorithmically. However, the simplicity of the "beads and springs" approach means that it is often difficult to correctly capture torsional rigidities, or directional (as opposed to symmetric) hinges. This becomes more of a limitation for ultra-coarse-grained models involving representations at lower resolution than the backbone C α atoms. An advantage of continuum mechanics models is that the material space between nodes is "filled-in" so that excluded volume is represented explicitly, and modeling torsional rigidity is simplified (assuming the appropriate torsional parameters can be obtained). However, in spite of the continuum approximation, there is still an inherent coupling between the discretization of the mesh, the material parameters, and the biomolecular flexibility in the model. While this provides the freedom to subtly optimize the system parameters to reproduce experimental results, it also provides an additional source of variability that needs a detailed physical understanding before an automated parameterization of such continuum models can be used reliably. Currently, considerable user intervention is required for each protein modeled, at every length scale.

| Nucleotide driven conformational changes and the chemo-mechanical cycle
Cytoskeletal molecular motors change conformation due to their interactions with ATP, their cargo, and the track. It is the precise choreography of these interactions that gives rise to the chemo-mechanical cycle underlying the biological function of the motor, including directionality. In atomistic simulations, there is sufficient detail that conformational change can (in principle) be predicted from the simulations, but in practice these timescales are orders of magnitude too long for even a single step in the cycle to be captured.
If atoms are not represented, then the nucleotide is usually too small to be included explicitly and ATP/GTP driven structural changes need to be imposed by the modeler (as in Refs. 63,89). For many cytoskeletal motors, X-ray crystallography and cryoEM/ET have provided detailed information about some of the intermediate structures that the motors adopt during each step of their chemo-mechanical force cycles, but key information is missing. For example, in myosins transient states such as the ADP-P i bound pre-powerstroke state are hard to capture. Although there is a wealth of crystallography data using stable analogues such as ADP.VO 4 to mimic these, none are filament bound. 30,103,122,123 Filament bound states can only be captured via cryoEM and all current structures are in the post-powerstroke/Apo state, 28,30 as the pre-powerstroke state is short lived and has weaker binding affinity. The kinetic rates associated with each step have been measured using biochemical assays 118,119,124 and microscopy, 41,52 and so can be used to inform computer models. These rates theoretically depend upon the change in mechanical energy, and hence the motor configuration, between steps, 125 which implies that structural deformations due to thermal fluctuations may affect the probability that a conformational change will occur.

| Track-motor interactions
Changes in affinity are absolutely crucial in the action of molecular motors because they choreograph the chemomechanical cycle. For example, in dimeric Myo5a, both motors are bound to the track in the strained high affinity ADP-bound state. 24,126 Biased ADP release in the rear motor allows a new ATP to bind and the motor to detach, releasing the tension on the front head and enabling the rear head to move forwards. 24,118,119 An equivalent mechanism has been proposed for cytoplasmic dynein. 79 While biomolecular affinities, such as track-motor interactions can often be quantified experimentally, where these are not available it is necessary to estimate them. At the atomistic level, all interatomic energetic interactions can be included explicitly, but the entropic contribution is largely missing, because of limited sampling timescales. Therefore, even though atomistic models in principle contain all the necessary information to predict biomolecular binding affinities (assuming perfect MD force-fields), in practice this remains one of the big unmet challenges for computational modeling.
Free energy landscapes become "smoothed" at the mesoscale, because these entropic contributions are implicit in the interaction potentials. For coarse-grained and ultra-coarse-grained particle-based MD, specific intermolecular interactions are computed "bottom up" from a potential of mean force (see Ref. 127 for details of the theory and Ref. 128 for an accessible overview) using a more detailed (usually atomistic) MD simulation. This procedure is extremely computationally expensive, because the entropic component is very slow to converge, and the potentials are not easily transferrable to other situations. Alternatively, new computational methods for computing biomolecular entropies, that systematically consider a multiscale hierarchy of dynamical interactions, from bond vibrations through to whole protein translations and rotations, may provide sufficient theoretical understanding that we can provide estimates of these contributions from an empirical model. 129 The representation of biomolecular affinities at the mesoscale is highly challenging conceptually, because it depends on the level of detail. Compared to the atomistic level, different types of motion emerge at the mesoscale. The paradigm at the atomistic level is highly specific molecular recognition, because the dominant techniques (e.g., X-ray crystallography) for structure determination select the unique minimum energy conformation. Specific interactions and binding affinities can also be defined and incorporated based on experimental data. However, at the mesoscale, additional types of affinities and related dynamics emerge, and non-specific interactions play a larger role, as was observed using coarsegrained models to understand diffusion on posttranslationally modified MT surfaces (Section 3.1). 39,40 One example is molecular sliding or gliding, in which non-specific interactions reduce the diffusional space of the head of a motor from 3D to quasi-1D. 39 Rather than defining a specific binding site, here the essential physics that needs to be captured is a diffuse interaction that is "smooth" in either one or two dimensions, allowing free diffusion, but restricted in others. This requires fundamentally different interaction potentials at the mesoscale to atomistic models.

| Mesoscale biomolecular dynamics and the surrounding environment
Viscous fluid media are ubiquitous in biology, be it the cytoplasm, nucleoplasm, blood plasma, or synovia, within which most (if not all) biomolecules and organelles are suspended. Hence any realistic description of dynamics of biomolecules must account for the huge extent to which motion is slowed by viscous drag. All biological fluids are inhomogeneous, but when modeled at a sufficiently coarse-grained scale can be approximated as a homogeneous fluid. For example, at the small scale of 10 nm, the cytoplasm is a crowded and inhomogeneous suspension of proteins, macromolecular complexes and other biomolecules in cytosol, 96,130 but at a scale of 10 μm may be considered as a uniform viscous fluid, which greatly simplifies the problem. 131,132 The viscoelasticity of the cytoplasm, a combination of instantaneous (elastic) and time-dependent (viscous) responses to deformation, and molecular crowding are further considerations that can profoundly impact intracellular biomolecular dynamics, 133,134 but are not discussed in this review.
One of the simplest effects of viscosity is the heavy damping of motion of immersed objects due to viscous drag. The fluid dynamics of most mesoscale systems corresponds to the low Reynolds number, or Stokes' limit, meaning viscous forces, that is, those arising from intermolecular forces in the fluid, dominate over inertia. 135 We consider here the basic example of a spherical particle diffusing through homogeneous cytoplasm, subject to a viscous drag force. The diffusivity of the sphere is given by the Stokes-Einstein equation where k B is Boltzmann's constant and R is the hydrodynamic radius of the sphere. This gives a characteristic (average) diffusion time for traveling a linear distance L of So, comparing a small, globular protein such as ubiquitin (R Ub $ 2 nm 136 ) with a vesicle (R Ve $ 100 nm), for a temperature T = 310 K and fluid viscosity η = 10 À3 Pa s 137 we obtain diffusivities D Ub = 100 μm 2 s À1 and D Ve = 2 μm 2 s À1 respectively, meaning they diffuse a distance of 1 μm in around τ Ub = 5 ms and τ Ve = 250 ms. To confirm this inertialess description of thermal diffusion we can consider the timescale, τ I , over which the inertia of a particle decays (i.e., inertia has a negligible effect at timescales greater than this), given by where m is the particle mass and ρ $ 10 3 kg m À3 is the mass density of a typical protein. 138 This yields τ I,Ub = 10 À12 s and τ I,Ve = 10 À9 s, which when used in Equation (2) gives the length scales over which inertia decays as L I,Ub = 10 pm and L I,Ve = 60 pm, comparable to atomic bond lengths. 139 We can hence neglect the effects of macromolecular inertia for diffusion at the mesoscale.
If we now consider motion that is directed rather than random, as is the case for cytoskeletal molecular motors and cytoplasmic streaming flows, then the distance traveled can be approximated as L = v τ, where v is the speed and can be taken as $0.5 μm s À1 for cargo-laden motors in animal cells. 132 Substituting this into Equation (2) and rearranging, we get τ = 2D/v 2 , the time scale beyond which motor transport is faster than diffusion, which yields τ Ub = 800 s and τ Ve = 16 s; the corresponding expression for length scale, L = 2D/v, gives L Ub = 400 μm and L Ve = 8 μm. In general, random and directionless diffusion is sufficient when the protein radius or length and time scales are small (L / t 1/2 ), but size-dependent viscous drag reduces diffusivity and makes directed transport via motors and flows preferable for larger macromolecules (L / t). For example, in Lombardo et al. 4 groups of Myo5a carry a 350 nm vesicle at a speed of 0.4 μm s À1 . Based on the above, motor transport is faster than diffusion across distances greater than 2.8 μm, which is less than the distance recorded in the study (3.7 μm), 4 demonstrating motor transport is advantageous.
One broad and actively studied class of intracellular flows is cytoplasmic streaming: directed flows occurring in large eukaryotic cells such as Drosophila oocytes (hundreds of μm in diameter), with flow speeds ranging from 25 to 300 nm s À1 , 140 and Chara algae internodes (up to 10 cm in length) with maximum speeds of 100 μm s À1 . 141 These flows, which vary widely in magnitude and patterning across cell types, accomplish nonselective transport of groups of macromolecules and organelles across large expanses of a cell, acting as an additional method of directed transport alongside cytoskeletal motors. 132 Streaming is involved in cytoplasmic mixing 142 and phase segregation during meiosis 143 as well as the positioning of organelles, for example, pushing the meiotic spindle towards the oocyte center in mice 144 and correctly polarizing MT networks required for Kinesin-1 transport in Drosophila oocytes. 145 Despite occurring across potentially vast distances, streaming originates in part from the mesoscale activity of cytoskeletal motors: cortical actomyosin contraction governed by Myosin 2 and the Arp2/3 complex, 144 and hydrodynamic coupling of cargo carried by Myosin 11 (a plant analogue of Myosin 5 146 ) 147 and Kinesin-1 140 to the cytoplasm have all been shown to generate streaming flows.

| Modeling viscous effects
The drag force, F, on a macromolecule is given by F = ζv, where ζ is the drag coefficient that depends on its size and shape, but is approximately proportional to the longest length, l, and v is the velocity relative to the background flow velocity. For globular proteins and vesicles, coarse-grained models tend to use the spherical coefficient shown in Equation (1) for the sake of simplicity 63,86,90,140 ; high aspect ratio particles such as rods have an anisotropic drag where the resistance to motion parallel and perpendicular differs by up to a factor two, although some models ignore this and use a single coefficient. 82,84 Hydrodynamic interactions (HI) refer to the fluid-mediated, correlated motions that occur between macromolecules in viscous fluids. When a macromolecule moves it induces a flow velocity field within its local environment, the effective range of which decays as l/r, where r is the distance from the macromolecule. 148 Any other macromolecules within range will be perturbed by this flow, which in turn causes them to induce local flows of their own; HI can thus propagate through a biomolecular system and affect the overall dynamics in a nonlinear fashion. Studies have demonstrated that HI lead to synchronized motion in closely packed systems of filaments such as cilia arrays 149 and bacterial flagella bundles, 150 accelerate the folding rates of proteins, 151 and are critical to the generation of fast flows during cytoplasmic streaming. 140 The long range and complex nature of HI makes them computationally expensive to model: N-particle Brownian dynamics scales as O(N 3 ) when including HI because the dynamics of each particle and the stochastic noise due to Brownian motion are coupled to every other particle in the system. 148,152 Spatially resolved models that include steric interactions without explicit solvent may need to include HI as well as localized viscous drag. Where filaments sterically interact they are also likely fall within the l/r effective range of HI, unless the environment is sufficiently packed so as to screen the interaction and shorten its range. 153 Coarse grained models are necessary to reach the relatively slow timescales at which flows occur and to compensate for the algorithmic complexity of HI.
The use of Brownian dynamics to study both dry (no HI) and wet (with HI) active systems like the cytoskeleton or bacterial suspensions is widespread, along with many other well-established techniques that consider flow and HI at various levels of detail. [154][155][156] The stochastic immersed boundary method developed by Atzberger et al. 157 captures HI, thermal fluctuations, and structural deformations of mesoscale macromolecules. Simulations over 3 ms of a simplistic 1D motor transporting a meshed 3D vesicle under flow reproduced the drag force of a spherical vesicle, although the method is general and applicable to systems with greater degrees of freedom. The model for cytoplasmic streaming in Drosophila oocytes by Monteith et al. 140 treated MTs tethered to the cortex as hydrodynamically interacting chains of spheres, upon which sat additional spheres representing Kinesin-1-cargo complexes with adjustable separations that imparted forces onto their parent MT; the model reproduced the slow and fast streaming regimes observed by confocal microscopy over length and time scales of 100 μm and 800 s, respectively. The 2D whole-cell model of migration by Madzvamuse et al. 94 represented the actin network as a viscous gel with forces being imparted on the intracellular fluid from the contractility of Myosin 2 and polymerization of actin, both of which were represented implicitly by concentrations that moved via diffusion and advection in the absence of HI. The resulting polymerization, flow velocity, and cellular expansion in the direction of an actin gradient, along with an observed linear relationship between expansion and total amount of actin, highlighted actomyosin as a significant factor in driving cellular migration.

| VISUALIZATION AND MODEL BUILDING FOR THE BIOLOGICAL MESOSCALE
MD has benefited from visualization tools for presenting and understanding structural biology for decades. [158][159][160] However, as the field moves into the mesoscale, it is time to consider what tools will be needed for researchers to visualize their data.
Biomolecular simulations follow a workflow with visualization tools used to help with each of the steps involved. A system of connected particles or a mesh is built and then parameterized to create a model, a simulation using the model(s) is run (some models are used for static analysis) and the simulation results are visualized for analysis. Finally results may be visualized to present to the relevant audience such as in papers or to the public. This section discusses visualization tools designed for this workflow, many of which are integrated tools that can be applied to multiple steps, and have functionality suitable for studying the mesoscale.

| Visualization for model building
Particle and mesh-based models are two paradigms, not normally used together, which are needed to simulate and visualize the mesoscale (Figure 3). Both require detailed shape information in order to construct the coarse-grained model of the system. UCSF Chimera (Figure 4) is a molecular graphics system designed to analyze molecular structures with extensions including viewing 3D grid data such as electron density maps, viewing movies of trajectories, and displaying sequence alignments, and facility for users to develop their own extensions. 164 Chimera is commonly used to visualize cryoEM maps and rigid fit atomistic models to them, 165,166 therefore this could provide an interface for the particle-and meshbased paradigm. Recently the developers have released ChimeraX (Figure 5), which aims to improve performance and user experience; support for modern techniques in areas such as virtual reality (VR), light-sheet microscopy, and medical imaging data; and more tools for external developers such as an app store for extensions. The software is still in its infancy so there are multiple areas for future development, including further analytical and developer tools. 167 One tool that enables visualization and assembly of meshed macromolecules is GAMer, a tool that converts structural datasets from cryoEM or atomistic information (e.g., PDB files) into geometric meshes used in simulation and/or visualization. While primarily in C++, GAMer also has a Python application programming interface, as well as an add-on for Blender (Figure 3, right) so users have a GUI, allowing them to visualize and edit surface meshes. 163 It is clear that the functionality of this tool will be useful in the future of integrative modeling.
A bottleneck in the extraction of macromolecules from imaging methods such as cryoEM is the need to reconcile the high level of structural complexity of the molecules, with the limitations of imaging techniques. 169 The mesoscale exacerbates these issues because as sizes approach the super-macromolecular the level of complexity increases, especially when macromolecules become components of a whole organelle or cell. 158 An issue with building subcellular architectures is that with the component macromolecules it is difficult to know their quantities, orientation, and F I G U R E 3 (Left) The particle-based paradigm is represented by the atomistic structure of ADP-BeF 3 bound myosin 5a (PDB: 1W7J 161 ) zoomed into the active site in visual molecular dynamics (VMD). 162 (Right) The mesh paradigm is represented by a surface mesh of dimeric heavy meromyosin 5a (as seen in Figure 2d) generated using the FFEA toolset, 86 visualized in Blender (https://www.blender.org/), and edited and tetrahedralised using the BlendGAMer addon 163 (panel to the model's immediate right) F I G U R E 4 EM electron density map (gray) (EMDB: EMD-23082) with fitted model of axonemal dynein (orange) bound to a microtubule (blue) from 29 (PDB: 7KZM), visualized in UCSF Chimera 164 distribution. Tools are needed for analyzing experimental microscopy data to try and understand subcellular architectures better, to get a clearer picture of how these components are distributed, and to be able to share these results for collaboration in a system such as PDB-Dev (see Section 6).
While progress has been made in automating the setup of mesoscale simulations, there are other methods worth mentioning too. Most importantly, the increased use of machine vision for processing cryoEM data [169][170][171][172] introduces the possibility of integrating similar techniques into the visualization system. For example, creating a tool that takes EMDB data, identifies components in it (e.g., MTs and proteins) and passes them to the simulation. This could be further aided by the recent advancement of creating ENMs from EM maps to gain a series of conformers. 81

| Visualization of simulation results
For multiscale simulations the approximate shape of the molecule must be generated. [158][159][160] As discussed in Section 4.1, there is a lower size limit for objects before they are made visible. However, objects that are represented implicitly may also need a visual indication of their spatial distribution. Different background colors could be used to represent the chemical makeup of the solvent, for example where spatially resolved models are coupled to reaction-diffusion equations. 19 Often the visualizations of mesh-based simulations use pseudo-color to show the magnitude of the simulated properties. 94,173 For example, using a color gradient to show the magnitude of a force that is applied to a node, or using colormaps to represent different parameters such as temperature or material parameters. The use of color needs careful consideration to ensure that the visualizations are easy to understand and are accessible as choosing certain colors may obfuscate the results for people with defects in color vision. 174 It is also important to consider how the visualization system shows the system boundary conditions, for example so that users can detect artifacts due to finite size effects.
Flow is another important characteristic of simulations that would be useful to visualize. It is often desirable to see which parts of the protein were deformed or moved due to flow, and the speed and the correlation between flows in different regions of the fluid. A simple solution to this is to have the visualization system present the results as an F I G U R E 5 Atomistic model of an axonemal doublet MT built from single particle cryoEM (PDB: 6U42 168 ) visualized in UCSF ChimeraX 167 animation of the simulation at different time steps, but showing complex motion in static images such as those needed for papers is more difficult. MoFlow, for example, depicts flows using of path lines and ribbon surfaces. 175

| VR, stereoscopic projection, and interactivity
It could be easy to dismiss the use of immersive technologies in research as a gimmick, but for biomolecular simulations, which involve numerous 3D objects moving in space, it can be advantageous to explore and manipulate the system. Some have developed tools such as the Interactive Molecular Dynamics in Virtual Reality (iMD-VR) where, the focus is on education and science communication, allowing people to have a fully interactive stereoscopic (3D) view of a working system as opposed to more traditional methods of 2D textbook diagrams and pre-recorded video. 176,177 As for research, there have been cases of groups using VR for effective study, such as using iMD-VR and Narupa to interactively and flexibly dock drugs into binding pockets, 178 and using ensemble MD simulations and VR visualization (eMD-VR) to study ms timescale protein conformational changes. 179 It is easy to see how this could be useful for interactively building super-macromolecular arrays and visualizing complex subcellular dynamics.
An approach for interacting with mesoscale systems is to use "painting" software techniques. CellPAINT is an illustration tool for creating complex scenes incorporating multiple structures from the PDB. It uses familiar "brush" tools to draw scenes, and stamp-like tools to attach components, such as motors to filaments ( Figure 6). There is a VR implementation for creating 3D scenes, and CellPAINT-Exo for enables hand-drawings to be imported. 180 LifeBrush is another VR tool for creating mesoscale simulation models, similarly using tools like "brushes" for creating scenes but completely in stereoscopic VR, because it is difficult using a 2D environment and a 2D input device to sketch a 3D structure. LifeBrush demonstrations (such as the LifeBrush video demo listed under Further Reading) hint F I G U R E 6 Simplified scene of actin with myosin 5a motors bound made in CellPAINT 180 using the default actin available and imported PDB 1OE9 15 towards how tools like these could be used to paint in features of a tomogram and build a simulation starting configuration from a reference image. The main challenge for experimental systems like CellPAINT and LifeBrush is that new software using non-standard inputs for building parameterized models requires part of the mesoscale simulation workflow to be software developer-lead, that is, every new molecular agent and its interactomics (e.g., interaction networks) needs new software to be developed. Increasingly sophisticated 3D VR toolkits of molecules, structures, and phenomena, with the ability for users to place them into their model, 181 will emerge as new molecules and functionality are added.

| CONCLUSION
The complexity of both protein 3D structures and the molecular recognition processes that drive the formation and activity of super-macromolecular protein complexes at the atomistic level has been well articulated. As new biophysical tools allow us to explore the biological mesoscale, we observe equivalent organizational complexity that is of equal importance to biological function. Biomolecular structures involving arrays of molecular motors are no exception. New community simulation and visualization tools that are multiscale and which can interface with disparate types of experimental data are required to provide a holistic view of the cytoskeleton, and cellular processes in general.
Public databases for computational models that use standardized file types that can be read by other simulation or visualization tools, akin to an extension of databases like the PDB (http://www.wwpdb.org/) and EMDB (https://www. ebi.ac.uk/pdbe/emdb/), will be key for integrating data, multiscale simulation, and visualization. PDB-Dev (https:// pdb-dev.wwpdb.org/) is a prototype archive that allows researchers to deposit integrative or hybrid models. 182 PDB-Dev provides a dictionary for integrative and hybrid models that allows for deposition of multi-state models containing ensembles of conformations, and multiscale models containing hierarchical representations that may be atoms, coarsegrained beads, or 3D Gaussians. Provision for inclusion of experimental restraints (e.g., chemical cross-linking data) and for descriptions of complex workflows combining experiment and simulations are also available. 183 Integrative modeling also requires interdisciplinary collaboration and cohesive scientific communities, particularly when agreement on file formats and code interoperability is required, 184 which will be aided by the development of collaborative visualization tools. The cytoskeleton provides an ideal biological system for stimulating such community activities.

CONFLICT OF INTEREST
The authors have declared no conflicts of interest for this article.

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no new data were created or analyzed in this study.