## Abstract

Molecular motors are proteins that excessively increase the efficiency of subcellular transport processes. They allow for cell division, nutrient transport and even macroscopic muscle movement. In order to understand the effect of motors in large biopolymer networks, e.g. the cytoskeleton, we require a suitable model of a molecular motor. In this contribution, we present such a model based on a geometrically exact beam finite-element formulation. We discuss the numerical model of a non-processive motor such as myosin II, which interacts with actin filaments. Based on experimental data and inspired by the theoretical understanding offered by the power-stroke model and the swinging-cross-bridge model, we parametrize our numerical model in order to achieve the effect that a physiological motor has on its cargo. To this end, we introduce the mechanical and mathematical foundations of the model, then discuss its calibration, prove its usefulness by conducting finite-element simulations of actin–myosin motility assays and assess the influence of motors on the rheology of semi-flexible biopolymer networks.

## 1. Introduction

### (a) Motivation

A defining characteristic of life is autonomous movement driven by metabolic activity. Complex macroscopic motion—from facial expressions to a javelin throw—results from the chemical and mechanical interaction of vast numbers of molecular motors and protein filaments. Their interplay makes *things move* [1]. In this paper, we present a method capable of elucidating fundamental mechanisms in active mesoscale, intracellular structures such as the actin cytoskeleton of eukaryotic cells. Potential applications include the study of force exertion in stress fibres and filopodia or the base unit of the myofibril, the sarcomere. Previously, we introduced a very powerful Brownian dynamics finite-element method for the simulation of cross-linked biopolymer networks where filaments and linker molecules are modelled with geometrically exact, nonlinear *Simo–Reissner* beam finite elements [2,3]. We also discussed methodic extensions that decouple the chemical topology of the filaments from their mechanical discretization [4]. Thereby, we laid the foundation for the present contribution, in which a computational model of a non-processive molecular motor is motivated, discussed and numerically validated. We focus on modelling motors interacting with semi-flexible filaments such as myosins, which enable skeletal muscle activity and important intracellular functions within non-muscle cells. Our intention is to provide a microscopic instrument whose purpose is the exertion of forces within large biopolymer networks. We aspire to use this instrument to study the mechanical effects of motor activity on these networks.

In the context of finite-element-based biopolymer network simulations, only very few cases model molecular motors discretely, e.g. in order to explain the strain stiffening in actin filament (F-actin) networks [5]. In this case, however, only the (athermal) filaments are modelled with beam elements, whereas motors are merely represented by force dipoles. So far, the finite-element method has been predominantly employed to model materials on the tissue scale where filaments and motor activity are just accounted for by special constitutive laws (e.g. [6,7]). Alternatively, molecular motors and their cargo have been resolved on much smaller length and time scales [8]. To the authors’ knowledge, no model so far offers a coarse-grained, yet discrete mechanical description of a molecular motor based on the finite-element method. Therefore, the primary goal of this publication is to close this gap. We believe that such a model is a crucial modelling component for the efficient simulation of large mesoscale network structures formed by individual, potentially cross-linked filaments.

In §1, we briefly introduce and discuss the most important motor proteins. Subsequently, we reprise the underlying Brownian dynamics finite-element approach and its extensions, which are employed in this publication (§2). Section 3 discusses the numerical model of a molecular motor based on a Simo–Reissner beam formulation. More precisely, the motor activity of a single myosin head domain serves as the archetype of the model. In §4, we present numerical examples, which calibrate the model and present further examples: simulations of *in vitro* motility assays are compared with experimental results. Finally, two rheological simulations of biopolymer networks (with and without motors) demonstrate the effect of our model on the viscoelastic behaviour of biopolymer networks. In §5, we provide the conclusion and a brief outlook on future applications.

### (b) Molecular machines

Motor proteins such as myosin or kinesin drive an abundant variety of intracellular processes. Among other tasks, they enable muscle contraction [9], vesicular transport and endocytosis, mitosis and meiosis (see [10]), or mechanosensing [11]. Molecular motors amplify and accelerate biological processes. According to Howard [9], the *cyclic* activity of motors can be characterized by three distance measures. The *working distance* *δ*_{w} describes the displacement per step, the *path distance* *δ*_{p} denotes the distance between consecutive motor-binding sites along a filament, and the *distance per ATP* *δ*_{ATP} is the displacement per hydrolytic reaction. These quantities can readily be converted into meaningful parameters of the motor model. Two other important quantities are the *cycle time* *τ*_{c}=*τ*_{on}+*τ*_{off} and the *duty ratio* *r*=*τ*_{on}/*τ*_{c}, denoting the duration of the entire cycle and the temporal fraction which the motor spends attached to its filament, respectively. The attachment time is *τ*_{on}, and *τ*_{off} is the complementary time. The duty ratio can vary greatly. Table 1 shows that *processive* motors such as kinesin feature values of *r*∼0.5 and above, i.e. they stay attached to a filament between individual steps. By contrast, most myosins are *non-processive*, detaching after each hydrolytic event, and have a low duty ratio *r*∼0.01–0.1 [15].

#### (i) Processive motors

Processive motor proteins of the kinesin and the dynein superfamily are the main molecular motors associated with microtubules. They transport proteins, nutrients and metabolites, but also entire organelles. Kinesins predominantly move axially parallel or helically towards the (+)-end; most dyneins move in the opposite direction [9].

#### (ii) Non-processive motors

Myosin motors bind to F-actin. They are propelled by ATP hydrolysis, which produces adenosine diphosphate (ADP) and a severed phosphate. Most myosins are non-processive, and predominantly move towards the *barbed* or (+)-end of the filament. They exist in muscle cells but also in other cell types, e.g. in cells of stereocilia [16] or the retina [17]. From here on, the term *myosin* will imply type II, which serves as the archetype for the motor model. Figure 1 shows the structure of myosin II with its two heavy chains (heads) and four light chains of amino acids.

The enzymatic cycle of myosin is divided into states (see [10]), which feature conformational changes of the heads [18]. In the initial, short-lived *rigour configuration*, one of the two heavy chain head domains is attached to the filament (figure 1*a*). When ATP binds to the head, it detaches from the filament due to a small conformational change and can slide along the actin filament (figure 1*b*). In figure 1*c*, a large conformational change triggered by ATP hydrolysis displaces the head in the (+)-direction. The head then attaches to the next binding site (figure 1*d*) and displaces the filament translationally by *δ*_{w}≈5–15 nm [12] depending on the myosin type. Myosin hydrolyses one ATP per step, i.e. its working distance is *δ*_{w}=*δ*_{ATP}. Having reattached to the filament, the severed phosphate is set free, triggering another large conformational change, the *power stroke*. The stroke is not explicitly depicted. Rather, it is the transition between the states of figure 1*e*,*f*. A major part of the cycle is attributed to the slow *recovery stroke* [19], which recreates the rigour configuration at the new position (figure 1*e*). Various experiments have found typical motor forces in the range of 1–10 pN [12,20–22] (see the electronic supplementary material, table S6).

## 2. Brownian dynamics simulations of biopolymer networks

In this section, we aim to give a compact summary of our Brownian dynamics finite-element approach employed in this paper, which we comprehensively described in [2,3] and which enables the simulation of large biopolymer networks [23] and their rheology [24]. It provides the basis for the new motor model discussed in §3.

### (a) Description of the simulation framework

A biopolymer network commonly has *three* major constituents governing its physics: *filaments*, cross-linking molecules (*linkers*) and a *fluid* phase, which contains the aforementioned molecules. The fluid commonly occupies a large volume fraction in intracellular environments. It is characterized by low flow velocities and low Reynolds numbers. Often, as in the present case, it is sufficient to model only its effect on structural components rather than considering complex hydrodynamic interactions [25,26]. We represent the other constituents *explicitly* and confine them within a cuboid computational volume cell with periodic boundary conditions as depicted in figure 2*a*. Semi-flexible filaments (e.g. F-actin) are the second constituent—slender structures with a length of a few micrometres. We model them as one-dimensional mechanical continua, which are discretized with Simo–Reissner beam finite elements [27]. Linkers represent the third constituent. They are also modelled as (short and very stiff) beams, and take on three distinct chemical states (figure 2*b*): free, singly bound and doubly bound. Here, we will introduce a fourth species: molecular motors, which interact with filaments in a way that is similar to that of linkers. In the most general case, all modelled constituents are initially distributed randomly within the volume cell. Brownian motion and chemical interactions between the constituents drive *self-organization*, which leads to fascinating network architectures [23].

### (b) Model of a semi-flexible filament

Semi-flexible filaments such as F-actin feature complex molecular structures (e.g. helices). We coarse-grain these filaments and model them as one-dimensional *Cosserat continua* with three translational and three rotational degrees of freedom (d.f.). According to Cyron & Wall [2,25,26], beams are an adequate mechanical model for such filaments, which can be discretized with nonlinear, geometrically exact beam finite elements. On a time interval [0;*T*] and along a curve parameter *s*∈[0;*L*], the geometry of a filament of length *L* can be described by a pair of functions
*a*and
*b*They determine the positions of the material points of the filament in ** θ** (see [27]). The pseudo-vector translates into a 3×3 triad

**,**

*t***,**

*n***} and employing the skew-symmetric pseudo-vector matrix**

*b**a*). The tangent

**is defined to point towards the (+)-end of the filament. Solutions to (**

*t***,**

*x***) are calculated by evaluating the beam's equations of motion**

*θ**a*and

*b*with internal elastic contributions (.)

_{el}and external contributions such as stochastic excitations (.)

_{stoch}, viscous effects (.)

_{visc}and other external loads (.)

_{e}.

*q*_{s}denotes the elastic section force. All forces and moments are given as loads per unit length. On the microscopic scale, inertial effects are negligible. The prime operator denotes (.)′=∂(.)/∂

*s*. Elastic forces and moments

*a*and

*b*Time derivatives are marked by superscripted dots. Matrices

*a*and

*b*with mean values 〈.〉, Boltzmann constant

*k*

_{B}, Dirac functions

*δ*

_{ss*}:=

*δ*(

*s*−

*s**) and

*δ*

_{tt*}:=

*δ*(

*t*−

*t**) and a constant temperature

*T*. Additional external forces and moments (e.g. steric or electrostatic interactions) are included in our code but are not considered here.

### (c) Model extensions

Filaments must provide a certain density of chemical binding sites along their geometry for motor–filament bonds. In [3,28], these binding sites coincide with the nodes of the beam finite elements. For motors, the nodal separation would have to equal the path distance *h*_{f}=*δ*_{p}. The gap between the previously chosen reasonable discretization length *h*_{f}=125 nm [23,24] and path distances as small as *δ*_{p}≈5.5 nm requires additional methodic efforts. Otherwise, the chemically required resolution leads to an insensibly fine mechanical discretization. Examples in §4b,c feature binding site separations 5.6–8.3 times smaller than the referential value given above.

The methods presented in [4] enable the arbitrary placement of chemical binding sites along the discretized filament, allowing for mechanical joints along the interpolated geometry, hence decreasing computation time by up to an order of magnitude. The map
*s* and the local (element) line parameter *ξ*∈[−1;1] relates positions on the filament to their position within a finite element. Next, we equip each binding site with an orientation by means of a fixed rotation relative to the material triad that defines the cross-section orientation at position *s*_{b} of the filament. At position *ξ*_{b}=*ξ*(*s*_{b}), this binding site orientation is given by the orthonormal triad
** Λ**(

*ξ*

_{b}) are updated with each simulated time step. Figure 3

*a*shows that the

*L*

^{2}-norm of the pseudo-vector

*θ*

_{Δb}=∥

*θ*_{Δb}∥ describes the relative rotation about

*θ*_{Δb}/

*θ*

_{Δb}=

*t*_{b}with respect to the second base vector

**of the interpolated triad**

*n***(**

*Λ**ξ*

_{b}). The binding site orientation

*n*_{b}(

*ξ*

_{b}) depends on

*θ*

_{Δb}, which can be used to map orientation patterns on the filament centre line [4]. Lastly, we define conoidal reaction volumes with orientations

*n*_{b}as shown in figure 3

*b*. The opening angle

*φ*

_{b}of the cone delimits the perspective of the binding sites. Employing this model extension, we determine which motor can bind to a filament. Whether a chemical bond may be established is determined by basic vector calculus and probabilistic measures (see §3c). Finally, the filament model requires a motor model which permits the establishment of mechanical joints at an arbitrary location along a filament. We will only briefly introduce this model as it is described comprehensively in [4].

## 3. Numerical model of a non-processive molecular motor

In this section, we will discuss a kinematically abstracted model of a non-processive molecular motor, the core novelty of this paper. It enables the exploration of effects of motor activity in simulations of structures on the scale of biopolymer networks. In analogy to the original Huxley model [29], the motor knows only two base states: *on* and *off*. For a power stroke, it performs a rotatory motion, which resembles a swinging cross-bridge or a lever arm [18,30].

### (a) General remarks

A coarse-grained and geometrically simplified representation of the molecular geometry suffices to study motor-related phenomena on the scale of filaments and above. We do not resolve the complex structural changes of active motors, yet we preserve the required kinematic and mechanical complexity to study, for example, active biopolymer networks. We propose the model of a contractile rotating rod, which allows for the principal representation of myosin motor activity. Motors in their power-stroke phase will be discretized with beam finite elements based on the formulation of Jelenic & Crisfield [27] and the extensions presented in [4]. Otherwise, only their position and chemical status is tracked.

The mechanism of force exertion is still controversially discussed, e.g. *thermal ratchet* models (e.g. [31]) in opposition to models where the power stroke itself is a large conformational change (e.g. [32,33]). Mixed models exist as well [34]. In our model, conformational changes are governed by different reference configurations, between which one may switch according to the current chemical status of the motor. The reaction kinetics will be accounted for in analogy to [3,4]. For now, a load dependency of forces and step lengths [22] is not included in the present model.

### (b) Dynamics of free motors

In general, as solutes of the surrounding fluid, motors are subject to Brownian motion. They move freely through the computational domain (figure 2). The balance equations for both molecule types for the *i*th discrete time step at time *t*^{i} read
*ζ*_{m}=6*πηR* [9], *R* being the size of the particle. Stochastic forces acting on the linker read
** I** denotes the second-order identity tensor.

### (c) Model of the actomyosin bond

The modelled motor has two reactive sites, one at each of its two ends. One can be interpreted as the head domain of heavy mero-myosin (HMM), which binds to filaments; the other is either fixed to a substrate or also free to interact with filaments. Experiments, with which we compare the simulations of §4, often feature only HMM without the light mero-myosin tail part (figure 1).

The reaction kinetics are modelled by Poisson processes with the association rate constant *k*_{+} and the dissociation rate constant *k*_{−}, which can be found in the literature (e.g. [35]). We use rate constants which are adapted to our reaction model:
*N*_{A} and modelled reaction volume *V*_{r} (for details, see [3]). The corresponding chemical association and dissociation potentials are given by
*k*_{off} is affected by tension or compression of magnitude *F* acting on the bond, which is modelled by Bell's equation
*k*_{off} in the absence of force and the force-dependent off-rate *k*_{bell} [36]. The characteristic distances Δ*x* are bond specific. Despite more complex approaches (see [37]), this model provides a good approximation of force-dependent chemical reactions. Our force measure is the magnitude of the axial force |*F*_{l}| computed within the motor beam element. Taking into account the axial strain *ε*, the force *F*_{l}, which enters equation (3.6), is
*x* determines the nature of a chemical bond:
*k*_{bell} increases with load. The bond between actin and myosin is characterized by Δ*x*=−2.5±0.6 nm [38]. Eventually, the modified probability of bond dissociation is given by

### (d) Cyclic motor activity and motor–filament interaction

The enzymatic cycle has to be reflected in the mechanical model. ATP is not explicitly simulated but is taken into account by means of a Poisson process, which models its reaction potential. The subsequent two steps, the separation of the phosphate and its release, which triggers the power stroke, are modelled as one step in the wake of the attachment process. There exist two basic conformations: a long state *δ*_{+}, the motor enters the drag stroke, which covers a distance *δ*_{−}=*δ*_{w}−*δ*_{+} where the attached motor is dragged along with the filament. Typically lasting for ≤1 ms, the drag stroke features a significant increase in the dissociation rate constant *k*_{off}∼2000 s^{−1} [9]. This way, a rapid detachment is ensured. The recovery stroke, during which the motor returns to state *τ*_{off} that determines the reattachment time scale. Since a recovering motor does not notably affect the transportation of filaments, it is not represented as a beam element. The recovery rate constant is given by
*p*(*i*_{r}Δ*t*≥*τ*_{off})=1 with *r*=0.05 (table 1), we get *τ*_{off}=38 ms (*τ*=40 ms [9]).

### (e) Polarity criterion

The establishment and the maintenance of a chemical bond between motor and filament also depend on the polarity of a filament, which results from the filament's anisotropy with respect to its longitudinal dimension. In the numerical model, we determine the polarity of a filament by means of the orientation of the binding site triad *Λ*_{b} from equation (2.8). We make use of its first base vector *t*_{b} as shown in figure 4, which is the unit tangent of the filament at the binding site. Polarity lets myosin bind only towards the (+)-end of the filament. We simply define that *t*_{b} points towards that end. The polarity condition reads
*t*_{b}. Here, *x*_{b} marks the location of the binding site and *x*_{m} is the point of reference of the motor. In figure 4*a*, the motor at position *x*_{m} is shown to be fixed to a substrate. The figure depicts two binding sites at positions *x*_{b,1} and *x*_{b,2} with their respective tangents *t*_{b,1} and *t*_{b,2}. Only the site at *x*_{b,2} is available for the formation of a chemical bond (figure 4*b*) due to the condition given by equation (3.12), i.e. *d*_{mb}∥*t*_{b}.

### (f) Geometric criteria

Apart from structural polarity, the model is required to satisfy additional criteria, which are motivated by geometry (→ motor size) and kinematics (→ possible motor configurations). First, a certain proximity of the motor and the considered binding site ∥*d*_{mb}∥=∥*x*_{b}−*x*_{m}∥∈[*L*_{m}−Δ*L*_{m};*L*_{m}+Δ*L*_{m}] is required. *L*_{m} denotes the size of the motor; Δ*L*_{m} is a tolerance. Second, we need to inhibit the formation of geometrically possible but physiologically highly unlikely bonds. To this end, the working distance *δ*_{w} helps to estimate an angular restriction of both the motor's and the binding site's reaction volumes (figure 4*b*). The method we use to create oriented reactive volumes [4] is employed to enforce this restriction. The size of the motor must confer with the specifications of myosin, which has a structure of length *d*_{l}≈8.5 nm serving the purpose of a lever [12]. During the power stroke, this lever rotates, and thereby covers a distance of approximately 10 nm [9]. The length of the entire head domain can be roughly approximated by 2–3 *d*_{l}. It provides an estimate for the motor size *L*_{m}.

The model is parametrized such that the working distance *δ*_{w} is achieved. In order to realize this working distance by contracting and rotating the motor beam element, the angle covered by the rod during the power stroke—from state *d*_{mb} and *d*_{mb,⊥} are the difference vector between motor position *x*_{m} and binding site location *x*_{b} and its projection orthogonal to the binding site tangent *t*_{b}, respectively. We additionally define a lower bound *φ*^{+}_{δ}/10). Allowing motors to bind to filaments within the interval *b*. The choice of these angles influences reaction kinetics as the reactive volume depends on them: making the interval too narrow diminishes chemical activity.

### (g) Modelling of conformational change

Force generation is triggered by the reaction model from §3d and is accomplished by an update of the translational and orientational reference configuration of the motor beam element. This update can be motivated in a very fundamental way. Conformational change is the effect observed during the translation from one local free energy minimum to another. The two states are separated by an energy barrier of some reaction-specific height. The energy needed to overcome that barrier is provided by ATP (or by thermal fluctuations depending on the model). Having overcome the barrier, the new conformation seeks to stabilize at its local free energy minimum. A change in the reference configuration of the beam model achieves an analogous behaviour. Equation (3.10) models the attachment and equation (3.10) provides the trigger for the change to another reference configuration, creating an off-balance restoring force that drives the beam to approach its new unstrained configuration, thereby displacing the filament. Hence, (re)setting the reference state not only has the same effect on the filament as a conformational change, but also adheres to the same energy principle.

The power stroke is modelled by a contraction *λ* and a rotation by an angle *t*_{b} at the motor attachment site. Such a modelling decision may differ from the more complex reality but entails a facilitation of the parametrization of the power stroke. Starting from state *a*), the configuration of state *t*_{b}. Our model generates an instantaneous contractile force by reducing the reference length *λ*=*d*_{mb,⊥}/*d*_{mb}. By a sensible choice of the binding angle interval [*φ*^{−}_{δ};*φ*^{+}_{δ}], the stroke distance *δ*_{+} of our motor model can be tuned to the experimentally observed stroke distance but may vary slightly depending on the geometrical situation.

Simultaneously, material triads *Λ*_{m} and *Λ*_{b} of the current configuration are used to compute the updated referential rotations of the beam element representing the motor. In figure 5*a*, the angle by which the triad needs to be rotated is given by
*d*_{mb} needs to be rotated to come to lie in the direction of *d*_{mb,⊥}. The rotation can be parametrized by its pseudo-vector representation

This modelling step reflects the idea of the swinging cross-bridge (see [18]). The applied beam formulation allows for the transduction of moments, which is why the incorporation of an updated reference rotation is of advantage. Otherwise, the longitudinal contraction has to work against the bending stiffness of the two joints. The entire update procedure is depicted in figure 5, which illustrates the mechanism leading to a translation of the filament. As a consequence of changing the reference orientation, the beam experiences restoring moments in *x*_{m} and *x*_{b}. Eventually, the superposition of contraction and rotation yields the stroke distance *δ*_{+} as illustrated in figure 5*b*. Using the inverse relation of equation (3.17), the update can be reversed if necessary, e.g. for the return to state *p*_{off} after the power stroke (see §3d).

The entire power-stroke sequence lasts for about 1 ms [9], of which a portion of *τ*_{−}≈0.4 ms can be attributed to the drag stroke. The remaining time *τ*_{+} represents the temporal fraction during which the motor is active. This time scale is used to quantify the rotation increment for a given time-step size Δ*t*, which simply is

## 4. Numerical examples

We present four numerical showcases, of which the first serves the purpose of motor calibration (§4a). The second example reproduces experimental results from [39] (§4b). The third example simulates a two-dimensional motility assay (§4c). The last example features a rheological comparison of two networks, one with and one without motors (§4d). Excluded volume or contact interactions between filaments are neglected here but will be included later based on a novel method modelling contact interactions for highly slender structures [40].

### (a) Motor calibration

First, we calibrate the motor in order to exert forces in the physiological range. A set-up as shown in figure 6*a* is used to assess the forces exerted by single-motor proteins [9]. A single-motor protein pulls on a filament, which is attached to a transducer element (e.g. a glass fibre).

The numerical example as shown in figure 6*b* mimics this set-up in order to record the force exerted by the motor beam element of length *L*_{m}. Another beam element acts as the filament to which the motor is attached. After the power stroke, which only entails an axial elongation of the filament, the motor element takes on a new state of mechanical equilibrium where we measure the axial force exerted onto the filament. We set the dissociation rate constant to *k*_{off}=0 s^{−1} and refer to the electronic supplementary material, table S1, for further simulation parameters. Except for its length *L*_{f}=0.1 μm, the properties of the filament element are listed in the electronic supplementary material, table S3.

Figure 6*c*,*d* shows the measured motor forces with respect to the two most important motor parameters, the stroke distance for the unloaded case *δ*_{+,0} and the size *L*_{m}. Exerted forces lie in the range of *F*∈[1;8] pN for parametrizations with *δ*_{+,0}∈[1;10] nm and *L*_{m}∈[5;25] nm, which are close to experimental results for myosin (see §1b(ii)). Due to the motor kinematics involving both translation and rotation, the increase in the exerted force with increasing working distance is nonlinear, as figure 6*c* shows. Figure 6*d* depicts the development of the motor force with respect to the motor size. Both figures provide data for fitting the motor model to experimental data. Depending on the motor species, its dimensions or force measurements, a specific parametrization of the motor beam element can be carried out.

### (b) Filament transport velocity in one-dimensional *in silico* motility assays

In order to validate our model, we chose a motility assay, a popular *in vitro* set-up, to study motor proteins under well-controllable conditions (e.g. [38,41,42]). Such an assay comprises motors and filaments, and can provide data on the transport velocity of filaments. In experiments, an assay usually consists of a micro-chamber with a volume of some tens of microlitres, which is filled with a fluid containing ATP. The motors are attached to a microscope coverslip. We compared our simulations of an assay with an experiment in [39] (plotted in [9], p. 221, fig. 13.3), which evaluates actin filament transport velocities *v*_{f} depending on the surface density *ρ*_{m} of myosin. A drastically increased viscosity (*η*=0.05–0.1 *Pas*) of the surrounding aqueous medium has no effect on the motor-induced filament transport velocity but diminishes lateral diffusion, thus facilitating the measurement of the filament velocity.

We assessed the *ρ*_{m}-dependent filament transport velocity of filaments in a one-dimensional problem. Three distinct configurations are depicted in figure 7*a*, in which a filament is transported in the direction of its (−)-end. As in the experiment, the filament's length was set to *L*_{f}=2 μm. It was discretized with 10 beam finite elements and featured binding sites every *d*_{b}=36 nm corresponding to the path distance *δ*_{p} of myosin. Here and in the following example, we calculated the filament velocity by averaging the magnitude of the velocities of the finite-element nodes over the course of the entire simulation. The motor surface density *ρ*_{m} translates to a corresponding line density ^{−2}. A periodic boundary condition was employed in the direction of filament translation at a repeat of *H*=10 μm. The temperature was set to *T*=303.15 *K* and the viscosity of the surrounding fluid to *η*=0.05 Pas [39]. The simulated time was *T*_{sim}=20 s at a step size of Δ*t*=5×10^{−4} s. The motor size was chosen to be *L*_{m}=20±10 nm, i.e. close to the size of a myosin head of 20–30 nm assuming a lever length of 8.5 nm [9] and a head of approximately double that size. The stroke distance in the unloaded case was *δ*_{+,0}=10 nm. The motor's cycle time was set to *τ*_{c}=40 ms with a duty ratio *r*=0.05. On the basis of experimental data [43], we set the motor model's circular cross-section to *A*_{m}=4.75×10^{−6} μm^{2} and Young's modulus to *E*_{m}=3.0×10^{6} pN μm^{−2}. A comprehensive list of parameters is given in the electronic supplementary material, tables S1 and S2.

With this parametrization, we achieved a good match to reported experimental values as documented by figure 7*b*. The *ρ*_{m}-dependent increase of *v*_{f} is captured well by the numerical model as well as the onset of saturation for *ρ*_{m}>1000 μm^{−2}. Mean values at high motor densities are 10–15% lower than the reported data points of the experiment. The difference may result from the simulation being strictly one dimensional, i.e. a real *in vitro* filament is displaced by more motors at a time than just by such lying exactly in the filament's trajectory.

### (c) Simulation of a two-dimensional motility assay

In this example, we employ the motor model to compare results with experimental values and to test the numerical feasibility of a dense two-dimensional field of motors transporting a multitude of filaments. Numerical comparisons with other state-of-the-art finite-element models are not possible as we present the first approach of this kind. The archetype of this example is an actin-HMM gliding assay, whose *in silico* realization is depicted in figure 8*a*. The filament density *ρ*_{f} characterizes an assay as either a *low-density* or a *high-density* assay. While filaments are disordered at low densities, they form collectively moving patterns at densities above a threshold *ρ**_{f}=20 μm^{−2} [44]. We examine and validate our model for the case of a low-density assay.

Since the experiment is two dimensional in principle, the simulated filaments are confined within the *xy*-plane. The visualization in figure 8*a* shows *N*_{f}=100 single filaments of length *L*_{f}=4 μm. The simulated area measures 6×6 μm^{2} and has periodic boundary conditions. Hence, the filament density is *ρ*_{f}≈3 μm^{−2}. Filaments are discretized with *N*_{e}=16 beam finite elements each with binding sites every *d*_{b}=*δ*_{p}=36 nm, and are initially straight. The temperature is set to *T*=293.15 *K*. The square matrix of motors features a distance *d*_{b} between two neighbouring motors, which translates to a motor density *ρ*_{m}=770 μm^{−2} (*N*_{m}=27 777 motors). The microscope coverslip is represented by an array of *substrate* filaments of length *L*_{s}=6 μm, to which the motors are attached, and each of which is discretized with *N*_{e,s}=20 elements, immobilized by Dirichlet boundary conditions. The motors are confined and spatially fixed in a plane parallel to the filaments’ plane of motion with an offset of 20 nm equal to the motor size *L*_{m}. The simulation covers a time interval *T*_{sim}=20 s at a step size of Δ*t*=5×10^{−4} s. Further simulation parameters, e.g. the filament's material properties, are compiled in the electronic supplementary material, tables S1 and S3.

We measure the transport velocity of filaments in analogy to §4b. Figure 8*c* depicts the distribution of filament velocities *v*_{f} with and without motors. Velocities were evaluated for time intervals of Δ*t*_{eval}=20Δ*t* for smoothing the motion of the filaments. A peak in the distribution can be seen at *v*_{f}≈3.5 μm s^{−1} followed by a smooth drawn-out descent towards higher velocities. The mean transport velocity is ^{−1}; see the electronic supplementary material, table S6, for literature values). As a reference, we evaluated the same filament ensemble in the absence of motors. The curve representing the case without motors in figure 8*c* shows a narrower distribution of filament velocities with the peak of the distribution located at approximately 1.6 μm s^{−1}. This is less than half the velocity found for the corresponding value of the motor-driven filament ensemble. The mean velocity is *in vitro* assay, the motor model provides directed acceleration through multiple motors cooperatively translating a filament.

Figure 8*d* shows the effect of forces acting on the chemical bond, i.e. the dissociation rate constant. The distribution has two major peaks surrounding the minimum of the experimentally determined HMM dissociation rate constant *k*_{−}=0.09 s^{−1}. The left peak relates to tensed motors during their power stroke. Due to employing Bell's equation (3.6), the dissociation rate constant drops by approximately 50%, making bond breaking less likely and increasing the number of motors moving the filament. The second peak marks an increase of *k*_{bell} by approximately 43% and corresponds to motors which are compressed before their power stroke. Motors which are in the long state

### (d) Motor-induced stiffening of cross-linked, semi-flexible networks

The final numerical example compares the results of rheological simulations with and without motors in analogy to [24]. To the authors’ knowledge, this simulated comparison is the first of its kind on this scale and on a broad frequency interval using an explicit motor model. The simulated network structure comprised 360 filaments (4 μm long, 16 elements per filament, binding site spacing *d*_{b}=62.5 nm), which were confined within a cubic, periodic volume with edges 6 μm long. The *N*=9000 additional non-filament molecules in the volume were divided into very stiff passive linkers and motors. The size of both passive linkers and motors was set to 100±20 nm (comparable to the cross-linking protein filamin used in experiments on active network stiffening [45]), an association rate constant *k*_{on}=90 s^{−1}, and a dissociation rate constant *k*_{off}=3 s^{−1}. The stroke distance of the motor was set to *δ*_{+,0}=20 nm. The reaction rate constant governing the power stroke was chosen in accordance with §4b. Apart from these motor-specific quantities, linkers and motors shared the same properties. All relevant parameters are listed in the electronic supplementary material, table S4.

After letting the network evolve for 100 s at a time step Δ*t*=2.5×10^{−3} s (figure 9*a*), we applied a unidirectional, sinusoidal shear at the top of the network with shear frequencies *f*∈[0.1 *Hz*,10^{4} *Hz*] (see the electronic supplementary material, table S5, for step sizes) and measured the stress response there. We fixed the bottom of the network by means of zero-displacement boundary conditions. The *same* network geometry was simulated with passive linkers (*N*_{l}=*N*=9000) and, separately, with *N*_{m}=900 motors and *N*_{l}=8100 passive linkers. Figure 9*b* depicts the rheological spectra gained by these simulations. Compared with the case with only passive linkers, 10% of the motors raise the storage modulus *G*′(*ω*) by up to 30%, thus increasing the frequency-dependent stiffness of the network as shown by the ratio of the storage moduli in the inset of figure 9*b*. Dissipation in the network increases by up to 75%, which is related to the increase in force-induced unbinding events due to motor activity. The amount of stiffening is not as excessive as reported in [45] due to the fact that we do not model myosin thick filaments comprising hundreds of myosins. It is comparable, however, to the stiffening of cross-linked networks with motors in the limit of rigid linkers in simulations with athermal filaments and force dipoles as motors (approx. 40% gain in stiffness, [5]).

## 5. Conclusion and outlook

In this paper, we have presented a numerical model of a non-processive molecular motor inspired by the motor protein myosin II, and which is based on a geometrically exact, nonlinear beam finite-element formulation. Our model accounts for conformational changes by means of changes in the reference configuration of the beam finite element representing the motor. Choosing a parametrization close to physiological conditions (motor stiffness, exerted force and stroke distance), our model offers a good approximation of single-motor behaviour and its interaction with filaments. It relies on very few model parameters and is therefore simple to fit to a multitude of non-processive motor species given that relevant properties are known, e.g. from experiments. The presented numerical model fulfils all requirements to be applied in complex studies involving the effect of motor proteins on subcellular network structures.

We have validated the motor model in a series of numerical experiments. We have calibrated the motor model in order to match the exerted force of the numerical model to its physiological counterpart. Subsequently, we tested the motor model in numerical representations of *in vitro* motility assays. We found close-to-experimental values for characteristic quantities such as the transport velocity of filaments in *in vitro* motility assays in one and two dimensions. To the authors’ knowledge, these exemplary simulations represent the first of their kind using beam finite elements. In simulations of cross-linked filament networks, our motor model induced internal stresses, stiffening these networks by up to 60% as compared with networks without motors. To the authors’ knowledge, we present the first simulation of the linear rheology of semi-flexible networks with motors on a broad frequency interval.

Future developments may include the extension of our model towards more complex structures, e.g. myosin thick filaments, which enable contraction of actin filament bundles [46–48] and alter the mechanics of actin networks [45,49]. Therein, a certain statistical dispersion in motor contraction rates, transport velocities and probably their sizes as well play a pronounced role. Such uncertainties may be included in the model, e.g. by assuming adequate probability distributions for these quantities. Another promising topic is the study of cooperative effects occurring due to motor activity, which have been observed in the contractile ring [50] and stress fibres [51]. Undoubtedly, robust and efficient numerical schemes for the contact interaction between filaments (and motors) is pivotal to a large range of potential applications in biophysics.

Finally, advanced formulations for mechanical joints between beams [52] can be used to emulate processive microtubule-related motor proteins such as kinesin and dynein. Such a step would further widen the scope, granting methodic access to other problems, e.g. the motor-driven intracellular transport of nutrients.

## Data accessibility

Simulations were conducted using the institute's multiphysics code. All presented data are reproducible employing the methods and results of our published work.

## Authors' contributions

All authors contributed to the design and the conception of the simulations. K.W.M. devised the computational method and conducted the simulations. K.W.M. and A.M.B. wrote the paper, which was critically revised by all authors

## Competing interests

We have no competing interests.

## Funding

This work was supported by the International Graduate School of Science and Engineering (IGSSE) of the Technische Universität München (TUM), Germany

## Acknowledgements

We thank Robijn F. Bruinsma and Alex J. Levine (UCLA) as well as Andreas Bausch and Oliver Lieleg (TUM) for fruitful discussions. We thank Christoph Meier and Maximilian Grill for helpful comments.

- Received August 12, 2015.
- Accepted November 24, 2015.

- © 2016 The Author(s)