F-actin bending facilitates net actomyosin contraction By inhibiting expansion with plus-end-located myosin motors

Contraction of actomyosin networks underpins important cellular processes including motility and division. The mechanical origin of actomyosin contraction is not fully-understood. We investigate whether contraction arises on the scale of individual filaments, without needing to invoke network-scale interactions. We derive discrete force-balance and continuum partial differential equations for two symmetric, semi-flexible actin filaments with an attached myosin motor. Assuming the system exists within a homogeneous background material, our method enables computation of the stress tensor, providing a measure of contractility. After deriving the model, we use a combination of asymptotic analysis and numerical solutions to show how F-actin bending facilitates contraction on the scale of two filaments. Rigid filaments exhibit polarity-reversal symmetry as the motor travels from the minus to plus-ends, such that contractile and expansive components cancel. Filament bending induces a geometric asymmetry that brings the filaments closer to parallel as a myosin motor approaches their plus-ends, decreasing the effective spring force opposing motor motion. The reduced spring force enables the motor to move faster close to filament plus-ends, which reduces expansive stress and gives rise to net contraction. Bending-induced geometric asymmetry provides both new understanding of actomyosin contraction mechanics, and a hypothesis that can be tested in experiments.


Introduction
The mechanics of actin filaments and myosin motor proteins in the cell cortex underpins movement (Yamada and Sixt 2019) and division (Pollard 2010) of biological cells. Early breakthroughs in understanding actomyosin dynamics occurred in muscle cells (Gautel 2011), in which actin and myosin form sarcomere structures. Sarcomeres involve filaments aligned in parallel with minus-ends in the centre and plus-ends pointing outwards. Relative motion of myosin motors towards filament plus-ends subsequently generates contraction by pulling filaments inwards. This mechanism is known as sliding filament theory (Huxley 2004). However, actomyosin networks in the cell cortex are disordered, with filaments distributed at random. Experiments (Murrell et al. 2015;Pollard and O'Shaughnessy 2019) and simulations (Tam et al. 2021;Ennomani et al. 2016) have shown that disordered actomyosin networks also contract (Chalut and Paluch 2016). According to sliding filament theory, filament pairs produce expansion if myosin motor proteins localise close to plus-ends, or contraction if myosin motor proteins localise close to minus-ends. In disordered actomyosin networks, motors localise near plus-ends and minus-ends with equal probability. Therefore, sliding filament theory alone cannot explain disordered network contraction. The origin of contraction in disordered actomyosin networks remains an active field of research.
Filament bending flexibility is commonly-hypothesised as a source of asymmetry that might explain contraction of disordered actomyosin networks (Murrell and Gardel 2012;De La Cruz and Gardel 2015;du Roure et al. 2019;Head et al. 2003;Tam et al. 2021). Actin filaments are semi-flexible (Stachowiak et al. 2014;Belmonte et al. 2017), such that they undergo small but significant bending (Broedersz and Mackintosh 2014;Murrell and Gardel 2012). Filament semi-flexibility is irrelevant in sarcomeres with parallel arrays of straight filaments, but is important for disordered networks in which motors can cross-link filaments at arbitrary angles and generate torque. Previous experimental and theoretical studies in disordered networks suggest that filament buckling gives rise to contraction. These studies show that filaments can sustain longitudinal tension, but buckle under longitudinal compression (Bidone et al. 2017;Belmonte et al. 2017;Cheffings et al. 2016;du Roure et al. 2019;Freedman et al. 2017Freedman et al. , 2018Lenz 2020;Murrell and Gardel 2012;Ronceray et al. 2016;Soares e Silva et al. 2011;Yu et al. 2018). This buckling mechanism can generate network-scale bias to contraction over expansion (Belmonte et al. 2017). Other studies have considered a related filament bending mechanism (Lenz 2014;Tam et al. 2021;Head et al. 2003;Popov et al. 2016;Kim 2015;Letort et al. 2015) as a source of force asymmetry. Bending involves applying forces that pluck filaments transversely, in contrast to the longitudinal forces involved with buckling. Lenz (2014) showed that filament bending produces forces that exceed those involved with longitudinal buckling, and Tam et al. (2021) showed that this mechanism facilitates network-scale contraction. A pertinent question is whether the force asymmetry provided by bending or buckling applies at the microscopic scale (Lenz 2014;Komianos and Papoian 2018;De La Cruz and Gardel 2015), or whether long-range effects transmit contractile force through the network, without requiring a microscopic asymmetry (Ronceray et al. 2016).
One approach to understand microscopic filament dynamics is to model a single filament as an inextensible elastic rod, as in 'worm-like chain' models (Broedersz and Mackintosh 2014;Lenz 2014). Broedersz and Mackintosh (Broedersz and Mackintosh 2014) used this approach to identify an asymmetry under extension and compression. Other authors have considered structures consisting of two-filaments and an attached motor (Lenz 2014;Belmonte et al. 2017;Hiraiwa and Salbreux 2016;Komianos and Papoian 2018). Lenz (2014) reported that disordered networks of rigid filaments with polarity-reversal symmetry (i.e. any configuration of filaments is equally likely as the same configuration with minus and plus-ends reversed) generate zero net contraction. Lenz (2014) also showed that filament bending gives rise to contraction in a twofilament system, and is the dominant mechanism of contraction for experimentallyfeasible parameters. During motor-induced deformation, semi-flexible actin filaments evolve to minimise their bending energy. Therefore, we model semi-flexible filament evolution as a curve-straightening flow. Mathematically, curve-straightening refers to deformation of curves in R 2 by decreasing their total squared curvature. Curvestraightening problems have been investigated since the 1980s Singer 1984, 1987;Linnér 1989Linnér , 2003. Wen (1993Wen ( , 1995 used the indicatrix representation and L 2 -gradient flow of the squared curvature functional to derive a fourth-order, semilinear parabolic partial differential equation (PDE) for the evolution of the curve. Oelz (2011) extended this work to model an open curve. However, theoretical analysis of curve-straightening flows is mostly limited to single curves, rather than the twocurve representations necessary to model a pair of filaments.
We extend previous curve-straightening models to derive a coupled PDE system for two symmetric, semi-flexible filaments with a myosin motor attached at their intersection. After obtaining the governing equations, we describe how to obtain the stress tensor, assuming the two filaments reside in a homogeneous background material. We then use asymptotic analysis and numerical solutions to provide a detailed explanation of how filament bending facilitates contraction on the two-filament scale. Our analysis suggests a contraction mechanism based neither on filament buckling, nor intrinsic force asymmetry where bending generates contraction. Instead, filament semi-flexibility creates a geometric asymmetry that inhibits expansion. Rigid filaments exhibit polarity-reversal symmetry, whereby contraction associated with a minus-end-located motor balances with expansion associated with a plus-end-located motor. Allowing filaments to bend breaks this symmetry, and the filaments become closer to parallel as the motor approaches the plus-ends. This decreases the spring force through the motor, enabling the motor to move faster close to the plus-ends. Fast motor motion inhibits expansive stress, and gives rise to net contraction. Our analysis provides a new hypothesis for bending-induced actomyosin contraction, and shows how contraction can occur on the microscopic two-filament scale. Fig. 1 Schematic representation of two actin filaments with a myosin motor attached at their intersection. Filaments are the curves z 1 and z 2 , and arrow heads represent minus-ends. The myosin motor protein is represented by the blue dot

Mathematical model
We develop a mathematical model for a myosin motor attached to two overlapping actin filaments. We represent filaments as open curves in R 2 , and denote their positions by z i (s(t), t) = (x i (s(t), t), y i (s(t), t)), for i = 1, 2 (see Fig. 1). The variable t denotes time, and s ∈ [0, L i ] is the arc length parameter, where L i is the length of the i-th filament. Actin filaments are polarised, so we adopt the convention that s = 0 corresponds to the filament minus-end, and s = L i corresponds to the plus-end. Since non-muscle myosin thick filaments are short compared to actin filaments (Dasbiswas et al. 2018), we represent the myosin motor as a point object existing at the intersection between the two filaments. We track its position by introducing the variables m i (t) ∈ [0, L i ], such that s = m i is the position of the motor head attached to the i-th filament. We assume that no other proteins, for example cross-linkers, are present.

Discrete force-balance equations
We express the mathematical model for the filament and motor mechanics as a system of force-balance equations, (2.1) The first three terms in (2.1) describe the drag, bending, and longitudinal stretching forces respectively on actin filaments. The fourth term represents longitudinal stretching along the myosin motor, and the final term describes forces between filaments and motors. We represent bending and stretching forces as the variation of potential energy, where terms involving δ denote variations. We formulate the force-balance equations (2.1) as a minimisation problem. This involves constructing a time-discrete scalar functional that contains a contribution for each force term in (2.1),  E a,stretch , and E m,stretch are the potential energies associated with filament bending, filament stretching, and motor stretching respectively. The terms E a,drag and E m,a are pseudo-energy terms with variations that correspond to finite-difference approximations of F a,drag and F m,a , which cannot be interpreted as variations of potential energy. The first term, E a,drag , describes drag friction between filaments and a passive background medium. Drag acts uniformly along the filaments and opposes filament motion. The term to represent drag between filaments and the background medium is where λ a is the filament drag coefficient, t is the time step size, and the superscript n refers to the previous time step in the discrete formulation, z n i = z i (s, t − t). We model bending of semi-flexible actin filaments via the elastic potential energy where κ a is the flexural rigidity, and primes denote differentiation with respect to arc-length, s. We assume that κ a is constant, and the same for both filaments. We obtain a term for filament stretching by assuming that actin filaments are inextensible.
To model this, we ensure that |z i | = 1 at every point along the filaments using the penalisation term where δ a is an arbitrarily small parameter that enforces the inextensibility constraints. The remaining two terms in (2.2) describe how motors contribute to the mechanics. To model motor stretching, we introduce another penalising potential, where δ m is an arbitrarily small parameter that penalises deviation from the constraint z 1 (m 1 , t) = z 2 (m 2 , t) stating that motors are to remain point objects. The final term in (2.2) describes interactions between filaments and motors. We assume that motors obey a linear (affine) force-velocity relation (Alcazar et al. 2019), illustrated in Fig. 2. This law integrates the multitude of force contributions exerted by myosin heads which decorate the myosin thick filament and interact with actin filaments (Kull and Endow 2013). Subject to zero force (when the motor is unextended), motors move with speed V m . As the motor extends, the spring force through the motor increases. We assume that motor velocity varies linearly with the spring force through the motor. If the force through the motor exceeds F s , the stall force, then the motor velocity is zero. The Fig. 2 The linear force-velocity relationship for myosin motors bound to actin filaments. An unextended motor subject to zero force moves towards filament plus-ends with the free-moving velocity, V m . If the spring force through the motor exceeds the stall force F s , the motor does not move corresponding pseudo-energy term consists of a linear term, and a quadratic drag-like term for the velocity reduction caused by the force through the motor, With the definition of (2.3-2.7), minimising the functional (2.2) for fixed t provides a time-implicit numerical method to solve the force-balance equations (2.1) for the filament and motor positions.
where primes denote differentiation with respect to arc length, dots represent derivatives with respect to time, and δ(·) is the Dirac delta function (not to be confused with variation). Equations (2.10) are a system of continuum force-balance equations for the filament and motor positions. They are formulated in a formal limit where δ a and δ m are small, and the force coefficients 1/δ a and 1/δ m in the variations of the penalising potentials (2.5) and (2.6) are replaced by the Lagrange multipliers λ 1 , λ 2 , and μ. Note that the sign of z 1 − z 2 in (2.10) will be absorbed by μ. As a consequence, solutions satisfy the constraints The equations are subject to the boundary and initial conditions where the subscript I represents an initial quantity. A detailed derivation of (2.10) and (2.12) is provided in Appendix A.

Calculation of forces and stress
An objective of this work is to describe how the filament and motor motion governed by our model generates contractile and expansive forces. We assume the two-filamentmotor structure is immersed in a dense network of cross-linked filaments covering a rectangular domain. A pair of actin filaments can locally deform the background network in which it is immersed. However, we assume the background can only undergo uniform elongation and shearing. A scenario supporting this assumption is that the background network consists of numerous two-filament-motor assemblies, all with the same shape as the reference pair that we describe explicitly (Fig. 3B). In this scenario, deformations occur equally everywhere in the domain, and the background network remains homogeneous. If the background is homogeneous, we can associate the tension at the boundaries of the domain with the stress being generated by the reference pair of actin filaments (Fig. 3A).

Fig. 3
Schematic representation of a two-filament-motor system existing within a rectangular homogeneous background medium. (A) The filaments are the curves z 1 and z 2 , and arrow heads represent filament minus (pointed) ends. Myosin motor proteins are represented by blue dots, and initially appear at the intersection between the two filaments. (B) Visualisation of the two-filament-motor structure embedded in a dense background network consisting of numerous assemblies (green lines) with the same shape as the reference pair (black lines) To quantify contraction, we compute the stress tensor for a small rectangular region of background material that encloses the two-filament-motor structure. The adjacent sides of the rectangle are given by the vectors L x = (L x x , L xy ) T , and L y = (L yx , L yy ) as shown in Fig. 3A. We compute the vectors F x = (F x x , F xy ) and F y = (F yx , F yy ), also shown in Fig. 3A. These vectors are the force components acting on the domain boundaries that must be applied to prevent uniform elongation and shear deformations.
These forces sum the contributions of both filaments and the motor, and provide a measure of net contractility. To derive closed form expressions for these forces, we return to the discrete formulation. We obtain F x and F y by first adding extra terms to the functional (2.2), and defining In (2.13), we use the modified drag term (2.15) The matrix L represents the transition from the coordinate frame L n x , L n y at time n, to the new coordinate frame L x , L y , corresponding to a rectangle that has undergone uniform shearing and elongation. If we impose L x = L n x and L y = L n y , the vectors F x = − ∂ L x E and F y = − ∂ L y E represent Lagrange multipliers that enforce the constant domain size and shape constraints. We assume that any possible deformations of the rectangle are small, such that Cauchy stress theory applies. The two-dimensional state of stress in the domain is then given by the stress tensor, (2.16) The bulk stress, then provides a measure of the contraction or expansion generated by the two-filamentmotor system. By convention, negative σ indicates contraction, and positive σ indicates expansion. The quantity σ is invariant to domain rotations, and equal to the average of the eigenvalues of σ . The associated eigenvectors of σ are the principal stress directions, which indicate the directions of maximum contraction or expansion.
To obtain an explicit expression for the bulk stress, σ, in terms of the filament positions, we differentiate the functional (2.13) with respect to L x and L y . This yields (2.18b) Applying the formal continuum limit t → 0, L → I, and (z i − z n i )/ t →ż i , we obtain Evaluating the bulk stress (2.17) then yields (2.20) Furthermore, the expressions (2.19) confirm that where z ⊥ i denotes a vector orthogonal to z i , and we obtain the result by substituting (2.10) forż i . The stress tensor (2.16) is thus symmetric, and the bulk stress σ is equal to the average of the eigenvalues of σ .

Nondimensionalisation
We nondimensionalise the PDE model (2.10)-(2.12) by introducing the length and time scalest = F s λ a L 2 a t, and x where hats represent dimensionless variables, and L a is a characteristic filament length. The dimensionless model is then (dropping hats for convenience) subject to the boundary and initial conditions and the constraints is a scaled Dirac delta function. The dimensionless parameters and forces are (2.26)

Model simplification
Before obtaining asymptotic and numerical results, we consider a simplification to the dimensionless model (2.23)-(2.25). First, we assume that the two filaments are symmetric about the vertical, that is and have identical length L a = L 1 = L 2 . Symmetry also implies that the relative position of the motor is the same for both filaments, m 1 = m 2 = m, and that λ * 1 = λ * 2 = λ * . To simplify the motor dynamics (2.23c)-(2.23d), we impose V m → ∞. Finally, we rewrite dimensionless flexural rigidity according to κ * = 1/ε, indicating that the flexural rigidity is large (ε 1), and that the filaments undergo only minor bending. On applying these simplifications, the dimensionless model (2.23) becomes (dropping asterisks on dimensionless parameters) subject to the boundary and initial conditions To obtain a measure of net stress, we integrate σ over the time between motor attachment and detachment. This yields The quantity J (T ) − J (0) describes the net, time-aggregated stress that the two filaments produce between motor attachment and detachment. This quantity will be important in our asymptotic and numerical investigation on how filament bending affects contraction.

Results and discussion
We analyse the simplified model derived in Sect. 2.5 to quantify how filament flexibility gives rise to contractile stress. First, we use asymptotic analysis to obtain an explicit approximation to the solution in the limit of infinite flexural rigidity. Through the leading-order problem, we show that a rigid two-filament-motor structure with polarity-reversal symmetry produces zero net stress. The first-order problem gives rise to a system of differential equations that governs the dynamics with small filament bending. Second, we obtain numerical solutions to validate the asymptotic analysis, and investigate solutions beyond the large flexural rigidity limit. These solutions reveal that contraction arises from a geometric asymmetry, whereby filaments become more parallel as the motor approaches the plus-ends. This inhibits expansion associated with plus-end-located myosin motors. Since contraction associated with minus-end-located motors is unaffected, the net outcome is a contractile two-filament-motor structure.

Asymptotic analysis
We construct an asymptotic approximation to the solution of the simplified symmetric model (2.28)-(2.30). Asymptotic analysis involves expanding variables in powers of ε, as ε → 0. On substituting the asymptotic series (3.1) into the model (2.28)-(2.30), the leading-order solution is the evolution of two rigid filaments with infinite resistance to bending. The first-order problem describes how small, non-zero bending affects the dynamics and stress. We present the key results and arguments in subsequent subsections, and give full details of the computations in Appendix B.

Leading-order solution
The leading-order solution describes the evolution of rigid filaments. To solve for z 0 , we consider the balance at O(1/ε). This yields The solution to the O(1/ε) problem (3.2) is a straight filament, whose direction we parameterise by the filament angle, θ/2, measured from the positive vertical axis (see Fig. 4). We write For a filament orthogonal to z 0 we use the notation where the symbol ⊥ denotes rotation to the left by π/2. A suitable ansatz for the position of a rigid filament solution satisfying the constraint (2.30b) is then where the leading-order motor relative position, m 0 , and leading-order vertical position of the intersection, y 0 , complete the parameterisation. The leading-order ansatz is illustrated in Fig. 4.
To obtain the leading-order solution, we consider the O(1) problem Schematic of a two-filament-motor system with rigid symmetric actin filaments. The myosin motor has relative position m 0 (t), and physical position (0, y 0 (t)). Filaments are symmetric about the dashed vertical line, which is the positive y-axis. The angle between the filaments is θ, such that the angle between a filament and the y-axis is θ/2.
Equation (3.6a)-(3.6b) are the leading-order equations governing the filaments and motor respectively. Equation (3.6c) provides two boundary conditions, and the solution must also satisfy the ansatz (3.6d) and orthogonality constraint (3.6e). To proceed, we use the orthogonality condition (3.6e) to infer the ansatz where h(s, t) is an arbitrary scalar function. Substituting the ansatzes (3.5) and (3.7) into the PDE for filament evolution (3.6a) enables us to solve for the leading-order quantities where H is the Heaviside step function, and ν 0 = m 0 − 1/2. Filament evolution then satisfies the ordinary differential equations where S = sin 2 (θ/2). Since z 0 is written in terms of the angle θ only, the system (3.9) determines z 0 . With h known, we subsequently obtain z 1 . Full details on this calculation are available in Appendix B.1. An important property of the system (3.9) is that it is invariant under a change of variables that reverses the direction of time. If we introduce the reversed-timet = T −t for an arbitrary constant T , we haveν 0 (t) = −ν 0 (T −t) andθ(t) = θ(T −t), and subsequentlyS = S. Consequently, if the motor is initially positioned at the pointed ends (ν 0 (0) = −1/2), and T denotes the time it reaches the barbed ends, then the time-aggregated stress vanishes, (3.10) This is because the equations and initial conditions satisfied by ν 0 , θ, andν 0 ,θ both coincide, and we have that ν 0 (t) =ν 0 (t) for allt ∈ [0, T ] (see also numerical result shown in Fig. 5C). This finding agrees with the previously reported results that rigid filaments with polarity-reversal symmetry produce zero net stress (Dasanayake et al. 2011;Lenz 2014).

First-order correction
The higher-order correction terms, z 1 , σ 1 and J 1 , elucidate the effect of small, nonzero bending on filament evolution and stress. To obtain expressions for the first-order corrections to bulk stress and J , we substitute the asymptotic expansions (3.1) into the stress terms (2.31) and (2.32). This yields Matching coefficients of ε then yields σ 1 = −2 1 0 z 1 2 + λ 1 ds, J 1 = 2 1 0 z 0 · z 1 ds. (3.12) In addition, we use the PDE (3.6a) and the ansatz (3.7) to obtain an explicit expression for the curvature of z 1 , Since the first-order correction to stress, σ 1 , involves the currently unknown λ 1 , progress requires consideration of the O(ε) problem, which is Obtaining the solution to (3.14) involves an intricate calculation based on the ansatz where A(t) is a (possibly time-dependent) constant of integration. The form of (3.15) arises from the ansatzes (3.14e) and (3.7), and gives rise to a system of equations for the degrees of freedom A(t), y 1 (t), and m 1 (t). We provide full details on the calculation to obtain this in Appendix B.2. A key result is the stress correction term, (3.17) Similar to the system (3.9), we can obtain a system of differential equations to solve for A(t), y 1 (t), and m 1 (t). Since h andh are in terms of the leading-order degrees of freedom θ and m 0 , we can subsequently compute σ 1 . However, the ODEs for A(t), y 1 (t), and m 1 (t) have no exact solution. Therefore, we continue our investigation using numerical solutions.

Numerical solutions
We compute numerical solutions to the simplified governing equations for filament and motor positions (2.28) in Julia. Our numerical method involves minimising the time-discrete functional (2.13) with t = 0.001. Energy minimisation is equivalent to a time-implicit numerical method for solving the dimensionless model (2.23). In our solutions, each filament has total length 1µm (Kamasaki et al. 2007), and consists of 50 equal-length line segments joined at nodes, about which segments can rotate. In Julia, we use the package Optim.jl (Mogensen and Risbeth 2018) to obtain the minimiser using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method. After obtaining the minimiser, at each time step we use automatic differentiation (ForwardDiff.jl) of the functional (2.13) to compute the forces F x and F y , and subsequently bulk stress σ (2.17).

Comparison with asymptotic analysis
We begin by computing numerical solutions for two symmetric filaments with m(0) = 0, and θ(0) = π/2. Like the asymptotic analysis, we assume these filaments are initially rigid, V m → ∞, and solve until the motor reaches the plus-end and detaches. First, we compute a solution for two rigid (ε = 1 × 10 −5 ) filaments, to validate the leading-order bulk stress σ 0 = 2ν 0 , and the solution to the system of ODEs (3.9), which governs z 0 . As Fig. 5A and B show, for both of these we obtain agreement between the numerical solution and leading-order solution. Furthermore, Fig. 5C illustrates the result from (3.10), namely that zero net stress is generated when a motor traverses two rigid filaments from the minus to plus-ends, i.e. J 0 (T ) = J 0 (0) = 0, where T = 0.627 is the time at which the motor reaches the plus end. Next, we solve the model with ε = 0.01 to validate the formulae for h and σ 1 , (3.13) and (3.16) respectively. The dynamics of the two filaments and motor are illustrated in Fig. 6. As part of the solution, we compute h using the asymptotic formula (3.13) and numerical values of θ and m, and compare this with the numerical value for the curvature, As Fig. 7 shows, we obtain agreement between the numerical and asymptotic results. The curvature formula (3.13) also reveals the shape that the two filaments adopt as they evolve (the qualitative pattern is easier to see in Fig. 10). Initially, the filaments adopt a convex shape, as the positive curvature in Fig. 7A shows. As the motor moves and pulls the filaments inwards, their shape changes to concave, as Figs. 7C-D show. When the motor approaches the plus-end, the filaments return to a convex shape. The asymptotic result for h remains accurate for up to ε ∼ O(1), before breaking down for ε ∼ O(10). We also use the numerical solution with ε = 0.01 to validate the formula for σ 1 , the first-order correction to bulk stress. At each time step, we compute the stress σ, and compare with the stress in a simulation with ε = 1 × 10 −4 , which we consider to be σ 0 for rigid filaments. We then approximate the first-order correction as σ 1 ≈ (σ − σ 0 )/ε, and present results in Fig. 8A. For most values of t, it holds that σ 1 > 0. In particular, larger positive values of σ 1 occur close to t = 0 and t = T , or m = 0 and m = 1. Fig. 8A is surprising, because it suggests the introduction of filament bending generates stresses that are biased to expansion. Similarly, as Fig. 8B shows, the quantity J 1 (T ) − J 1 (0) > 0, also suggesting net expansive bias. Based on this, one Numerical results for a solution with two rigid filaments (ε = 1 × 10 −5 ). A Comparison between the numerical bulk stress, and the leading-order approximation given by (3.8). B Numerical validation of the system of ODEs (3.9). C Numerical result for J (t), confirming that J (T ) − J (0) = 0 for rigid filaments might conclude that bending cannot facilitate microscopic-scale contraction. However, we have not yet accounted for the changes in filament geometry, and how they influence motor dynamics. Further simulations in Sect. 3.2.2 will reveal this more clearly, and confirm that bending does facilitate net microscopic-scale contraction.

Flexible filament solutions
We now consider numerical solutions beyond the ε 1 regime considered in the asymptotic analysis. These solutions are with the same conditions as Fig. 6, where the motor is initially at the minus-ends of two symmetric filaments. We then solve the model until the motor reaches the plus-ends. Results are presented in Fig. 9. The quantity J (t) measures the effect of ε on net stress, and is plotted in Fig. 9A. For rigid filaments, we showed that J (T ) − J (0) = 0, indicating zero net stress as the motor moved from the minus to the plus-ends. Since J (T ) decreases as ε increases, the introduction of filament bending facilitates bias to contraction. This contractile bias is despite the quantities σ 1 (T ) and J 1 (T ) being positive, as in Fig. 8. Indeed, Fig. 9B confirms that σ 1 > 0, with stress increasing with ε close to t = 0 and t = T .
Semi-flexible filaments facilitate net contraction because bending breaks the polarity-reversal symmetry, and the resulting geometry favours contraction. As Fig. 9C shows, with increasing ε, the myosin motor moves faster along the filaments and detaches earlier. The increase in motor speed is largest as the motor approaches the plus-ends, which Fig. 9B shows is associated with expansion. As the motor approaches the plus-ends, the semi-flexible filaments adopt a convex shape that brings them closer to parallel at their tips, as illustrated in Fig. 9D. This decreases the spring force through the motor, enabling it to move faster. Since the motor moves faster close to the plusends, the expansive component persists for shorter time than the contractile component. Consequently, the time-integrated stress J (T ) − J (0) decreases as ε increases.
The results in Fig. 9 are relevant for in vivo actin filaments, for which the parameters (Kamasaki et al. 2007;Gittes et al. 1993;Thoresen et al. 2011;Reichl et al. 2008;Oelz et al. 2015) estimated in Tam et al. (2021) give ε = 68.5. To further our analysis, we compute a numerical solution with ε = 68.5 and V m = 1, to investigate whether contraction persists after relaxing the assumption of infinite motor velocity. The evolution of these filaments is shown in Fig. 10. Despite the slower motor speed, the evolution qualitatively follows the prediction from Fig. 7. Filaments are initially convex, then become concave, and adopt a convex shape again as the motor approaches the plus-ends. As Fig. 10F shows, the two filaments are curved when the motor reaches the plus-ends and detaches. To rule out the possibility that relaxation to straight configuration produces expansion that cancels out net contraction, we continued the simulation after motor detachment, until the filaments were again straight. We plot σ (t) and J (t) in Fig. 11. Although relaxation to the straight configuration (shown in Fig. 11A) generates a small amount of expansive stress, Fig. 11C shows J (2) − J (0) < 0, suggesting net contraction. Thus, our proposed geometric mechanism for contraction remains relevant for realistic filament flexural rigidity and motor speed. Since actomyosin networks (for example those in the cortex) consist of many cross-linked two-filament assemblies, our mechanism provides a possible explanation for the microscopic origin of network-scale actomyosin contraction.

Conclusion
Understanding the origins of actomyosin contraction is an open problem in cellular biophysics, with implications for cell movement and division. In this paper, we presented a detailed investigation of how a two-filament-motor system generates microscopic contraction if the filaments are flexible. We first derived a partial differential equation model, and described a method of computing in-plane stress. We then applied asymptotic analysis to a symmetric system with infinite free-moving motor velocity. The leading-order solution showed that two rigid filaments do not generate net stress if the motor traverses the entire length of the filaments. However, the introduction of filament bending enables the two-filament-motor structure to generate net contraction. This is because bending breaks the polarity-reversal symmetry of rigid filaments. The resulting geometric asymmetry draws the plus-ends closer to parallel as the motor approaches. This facilitates faster motor movement when motors are close to filament plus-ends, and inhibits production of expansive stress. Our analysis confirms that the microscopic dynamics of symmetric filament pairs and motors can explain contraction. We expect that the same mechanism also favours contraction in non-symmetric filament-motor assemblies and that, consequently, macroscopic contraction in disordered networks could arise from the accumulation of multiple filament pairs, without the need for nonlinear amplification of contractile stress. Nevertheless, the question of how these results apply to disordered networks remains open. In disordered networks, filament pairs cross at arbitrary angle and position, and interact with an active background of other filaments, rather than the passive medium considered in this work. Tam et al. (2021) confirmed that disordered networks described by the biomechanical model for semi-flexible filaments and motors presented in this study do contract. Another potential approach modelling disordered network contraction is to derive a coarse-grained, continuum model based on the assumption of infinite filament length (Oelz 2014). This, and investigating how micro-scopic mechanics give rise to structures including stress fibres (Pellegrin and Mellor 2007) and the contractile ring (Kamasaki et al. 2007;Svitkina 2018), will be subjects of future work.

Appendix A: mathematical model derivation
In this Appendix, we present a detailed derivation of the PDE model (2.10) -(2.12). We derive the system of force-balance PDEs using the variational principle. The functional for the structure consisting of a myosin motor attached to two semi-flexible actin filaments is where z i (s, t) = (x i (s, t), y i (s, t)), for i = 1, 2, are the filament shapes and positions, m 1 (t), and m 2 (t) are the motor relative positions. According to the variational principle underlying the derivation, the solution at the next time step is such that the functional derivative with respect to all degrees of freedom vanishes, that is Evaluating each of (A2) gives rise to the variational equations For the motor evolution equations (A3c) and (A3d), we can immediately apply the continuum limit t → 0, for whichṁ i = (m i − m n i )/ t for i = 1, 2. In addition, we consider the limit δ m → 0 and replace the forces due to motor stretching by the force μ which represents the Lagrange multiplier for the constraint This yields the ordinary differential equations which are force-balance equations for the myosin motor positions m 1 (t) and m 2 (t).
The equations (A5) represent that unloaded motors (for which μ = 0) move at the free-moving velocity, V m . As the motor moves, it is exposed to stretching forces, with magnitude given by the Lagrange multiplier μ. The term involving the dot product is the projection of this force onto the direction of motor movement along the i-th filament. Assuming a linear force-velocity relation, the ratio of the term involving μ and the dot product to the stall force, F s , determines the reduction of motor speed due to stretching forces.
The differential equations (A11) and (A5), and conditions (A9) and (A10) as well as the constraints (A4) and (A6) then define a system of force-balance equations and boundary conditions that govern the evolution of two inextensible, semi-flexible filaments connected to an inextensible motor. This is the dimensional PDE model (2.10)-(2.12) given in the main text.

Appendix B: Asymptotic analysis
In this Appendix, we present the asymptotic analysis of Sect. 3.1 in more detail. We first outline the method used to solve the leading-order problem (3.6), and subsequently consider the O(ε) problem (3.14) for the first-order corrections.
We now derive the ordinary differential equations (3.9) for y 0 , θ, and ν 0 , the three degrees of freedom that govern the leading-order filament position, z 0 . On taking the time derivative and variation, the ansatz (3.5) implies that Integrating (3.6a) against δz 0 then yields Substituting (B4) into (B5) then gives Collecting the coefficients of δ y 0 , δθ and δm 0 , we obtain the system of differential equations (writing ν 0 = m 0 − 1/2) where csc(φ) = 1/ sin(φ). The equations (B7b) and (B7c) for θ and ν 0 are also independent of y 0 , suggesting that the solution is invariant to vertical translations. Furthermore, the trigonometric functions can be eliminated by writing S = sin 2 (θ/2), which yields This completes the derivation of equation (3.9).
In the case where z I ,0 = z I the term z I ,1 is minimal in L 2 . Then, the degrees of freedom m I ,1 , y I ,1 , and A I can be computed using 0 = 1 0 z I ,1 · δz I ,1 ds, where δz I ,1 = 0 δ y 1 − z I ,0 δm 1 + z ⊥ I ,0 (s − m I ,0 )δ A.