\documentclass[a4paper]{revtex4}  %showkeys, onecolumn, prl,showpacs,twocolumn, ,12pt
\usepackage{latexsym}
\usepackage[colorlinks=true,linkcolor=blue,anchorcolor=blue,citecolor=blue]{hyperref}
\usepackage[hmargin={1cm,1cm}, vmargin={2cm,2cm}]{geometry}
\usepackage{float}
\usepackage{graphicx,amsfonts,amsmath,amssymb,amsthm,amscd}
\usepackage{epsfig}
\usepackage{bm}
\usepackage{xr}
\usepackage{indentfirst}
%\usepackage{autobreak}
\usepackage{appendix}
\usepackage{changepage}
\usepackage{color}
%\usepackage{lineno}
%\usepackage{boxedminipage}






\begin{document}
\clearpage
\newpage
\pagebreak
%\vspace*{1cm}


%\begin{appendix}

\setcounter{figure}{0}
\setcounter{section}{0}
\setcounter{equation}{0}
\setcounter{table}{0}
\setcounter{page}{1}
%\renewcommand{\thepage}{A\arabic{page}}
\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\thesection}{S\arabic{section}}
\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thetable}{S\Roman{table}}
\begin{widetext}
	\maketitle
	\section*{\bf Supplemental Materials for \lq\lq Cooperation of myosin II in muscle contraction through nonlinear elasticity"}\label{appendex}
	\vspace*{0.3cm}
	\centerline{Beibei Shen and Yunxin Zhang}
	{\it Shanghai Key Laboratory for Contemporary Applied Mathematics, School of Mathematical Sciences, Fudan University, Shanghai 200433, China.}
	\vspace*{0.4cm}\\
	

\setcounter{subsection}{0}

{\centering\subsection {Nonlinear Elasticity of Single Myosin II}}

Experiments \cite{kaya2010nonlinear} measuring the elastic property of  myosin motors for a wide range of positive and negative strains demonstrate that myosin motors have nonlinear elasticity.
If motors are under positive strain, they can be described as \lq\lq stretched" and approximated as worm-like chains (WLC) \cite{fang2021biomechanical,fedosov2010systematic,li2005spectrin}. On the other hand, if motors are under negative strain, they act as \lq\lq drag" motors, and their stiffness is significantly lower compared to positively strained myosin motors.
Low stiffness minimizes drag of negatively strained myosin motors during muscle contraction at loaded conditions.
The nonlinear elasticity implies that active myosin motors with high stiffness is primarily responsible for the forward step generation.
It is crucial to consider the nonlinear nature of myosin elasticity in order to establish a connection between the internal structural changes of myosin and physiological force generation as well as filament sliding.      
The force-extension relationship of a bound motor $i$ depending on the strain, $y_{i}$, as follows:
\begin{equation}
	f(y_{i})=
	\left\{\begin{aligned}
		&\frac{k_{B} T}{L_{p}}\left[ \frac{1}{4(1-(y_{i}/L_{C}))^{2}}-\frac{1}{4}+\frac{y_{i}}{L_{C}}\right],&y_{i} > 0 ,\\
		&\frac{k_{\mathrm{m}}}{\gamma}[\exp (\gamma \cdot y_{i})-1],& y_{i}\le 0,\\
	\end{aligned}\right.
	\label{F_WLC_S}
\end{equation}
where $f(y_{i})$ is a continuous function with continuous derivatives, $k_{\mathrm{m}}=3k_{B}T/2L_{p}L_{C}$, $y_{i}$ represents the strain of a bound motor $i$, $k_{B}$ is Boltzmann's constant, $T$ is the absolute room temperature, and $\gamma$ is a parameter. $L_{C}, L_{p}$ represent contour length,  persistence length of the chain, respectively, $L_{C}>0, L_{p}>0$. $y_{i}/L_{C}$ represents extending ratio, and the range of it is $( 0,1)$. 
The strain energy of motor is dependent on the strain $y_{i}$, 
\begin{equation}
	E(y_{i})=\int_{0}^{y_{i}} f(s) d s=
	\left\{\begin{aligned}
		&\frac{k_{B} T L_{C}}{4 L_{p}} \left( \frac{y_{i}}{L_{C}}\right) ^2\left( \frac{1}{1-(y_{i}/L_{C})}+2\right) ,&y_{i} > 0 ,\\
		&\frac{k_{\mathrm{m}}}{\gamma}\left[ \frac{1}{\gamma}\exp (\gamma \cdot y_{i})-y_{i}-\frac{1}{\gamma}\right] ,& y_{i}\le 0.\\
	\end{aligned}\right.
	\label{E_WLC_S}
\end{equation}

\begin{figure}[htbp]
	\includegraphics[scale=0.5]{EF.eps}
	\caption{The nonlinear force and energy of a bound myosin II motor as functions of extension ratio $y_{i}/L_{c}$. \textbf{(a)} The force-extension relationship of the bound motor $i$ exhibits nonlinearity, which depends on the extension ratio $y_{i}/L_{c}$. \textbf{(b)} The energy of the bound myosin II motor $i$ varies with the extension ratio $y_{i}/L_{c}$.}
	\label{Forcestrain}
\end{figure}

{\centering\subsection {General features of the
myosin II motors model}}
\begin{figure}[htbp]
	\includegraphics[scale=0.155]{cartoon_of_multiple_motors1.eps}\\
	\caption{The {\it microstate} vector cycle is the sequence of transitions $\mathbf{n}^{\prime \prime} \rightarrow \mathbf{n}^{\prime \prime \prime} \rightarrow \mathbf{n}^{\prime \prime \prime \prime} \rightarrow \mathbf{n} \rightarrow \mathbf{n}^{\prime} \rightarrow \mathbf{n}^{\prime \prime}$.}
	\label{S1}
\end{figure}
We model a myosin half-filament that has $N_{0}$ total motors, each of which
can be inactive or active. The number of active motors is expressed as $N_{T}$ and depends on the load $F$, which can be represented as 
\begin{equation}
	N_{T}(F)=N_{0}\left[\left(1-a_{0}\right)+2 a_{0}\left(\frac{1}{1+e^{-b F / F_{0}}}-0.5\right)\right],
	\label{NT}
\end{equation}
where $N_{0}=294$ is the total number of active and inactive motors, and $F_{0}=480$ pN is the isometric tension. Parameters $a_{0}=0.95$ and $b=7.76$ are determined by fitting to activation data from Linari $et. al.$ \cite{piazzesi2007skeletal}. 

We consider both the
{\it microstate} and the {\it mesostate} of the system. The {\it microstate} is defined by the vector $\mathbf{n}=\left(n_{\rm pre}, n_{\rm mid}, n_{\rm post}, n_{\rm det}, n_{\rm rec}\right)$, where $n$ is the number of motors in the state given in the subscript. Each of these individual $N_{T}$ motors is undergoing its mechanochemical cycle and can be in any of five states. The {\it mesostate} $N$ is the number of bound motors, $N=n_{\rm pre}+n_{\rm mid}+n_{\rm post}$. 

To better describe the mechanism of myosin II motors, we define the {\it microstate} cycle as the following sequence of states:
\begin{equation}
	\begin{aligned}
		\mathbf{n} &=\left(n_{\rm pre}, n_{\rm  mid}, n_{\rm post}, n_{\rm det}, n_{\rm rec}\right), \\
		\mathbf{n}^{\prime} &=\left(n_{\rm pre}, n_{\rm mid}, n_{\rm post}-1, n_{\rm det}+1, n_{\rm rec}\right), \\
		\mathbf{n}^{\prime \prime} &=\left(n_{\rm pre}, n_{\rm mid}, n_{\rm post}-1, n_{\rm det}, n_{\rm rec}+1\right), \\
		\mathbf{n}^{\prime \prime \prime} &=\left(n_{\rm pre}+1, n_{\rm mid}, n_{\rm post}-1, n_{\rm det}, n_{\rm rec}\right), \\
		\mathbf{n}^{\prime \prime \prime \prime} &=\left(n_{\rm pre}, n_{\rm mid}+1, n_{\rm post}-1, n_{\rm det}, n_{\rm rec}\right) .
	\end{aligned}
\label{n5_S}
\end{equation}
The {\it microstate} cycle for {\bf n} is the sequence of transitions $ \mathbf{n}^{\prime \prime} \rightarrow \mathbf{n}^{\prime \prime \prime} \rightarrow \mathbf{n}^{\prime \prime \prime \prime} \rightarrow \mathbf{n} \rightarrow \mathbf{n}^{\prime}  \rightarrow \mathbf{n}^{\prime \prime}$. Across this cycle, different myosin II motors step through the individual transitions, one ATP molecule is hydrolyzed, and the actin filament is pulled a distance.

The variable ${\bf n}, {\bf n}'$, etc. each represent a {\it microstate} vector that can be any possible {\it microstate} of the system, thus serving as a dummy variable. But, within a single equation, once one of them is defined, the other four are automatically defined in relation to it according to Eq.~\eqref{n5_S}.

Earlier, it has been observed that several myosin motors must act cooperatively during such step generation.
Recent experimental results \cite{kaya2017coordinated} provide persuasive evidence that the stepwise  displacements of actin is mainly generated  by the coordination of power strokes among two myosin motors.
To increase the possibility of coordinated power stroke, several motors must be  briefly \lq\lq stay''  in either the {\it pre-power stroke} state accompanied by the first stage power stroke or the {\it mid-power stroke} state accompanied by the second stage power stroke. 
This further reveals that the molecular properties of  myosin motors are specifically tuned to perform cooperative force generation for efficient muscle contractions.

Based on the fundamental concept of the parallel cluster model \cite{erdmann2012stochastic}, we introduce a novel approach in which each {\it microstate} ${\bf n}''$ possesses a distinct set of strain values. This advancement enables us to capture the strain distribution among the bound motors, where the strain of a bound motor $i$ is
\begin{equation}
	y_{i}( {\bf n}'', F)=y_{i}^{0}({\bf n}'', F)+d_{i}({\bf n}''),
\label{yi_S}
\end{equation} 
where $d_{i}(\mathbf{n}'') =0,d^{\prime}$ and $d$ for a motor in the {\it pre-power stroke}, {\it mid-power stroke}, and {\it post-power stroke} states, respectively. 
The stroke is divided into two components, the second component is slightly smaller in reality, so $d^{\prime} \ge 4$ in general. \\

{\centering\subsection {Myosin II  motors transition in a first-in, first-out (FIFO) manner}}

This model is approximate in general, but accurate in the limit that myosin motors transition in a first-in-first-out (FIFO) manner. Here, we describe this limit.

Consider the sequence of {\it microstates} shown in Eq.~\eqref{n5_S} for any initial {\it microstate} ${{\bf n}''}$. We give each bound motor an index $i \in  \left\lbrace 1, . . . , N\right\rbrace $ that indicates the order in which each motor bound the actin filament, where $i = 1$ is the myosin that most recently bound to the {\it pre-power stroke} state. 

For FIFO transitions, the power stroke transition ${\bf {n}''' }\to{\bf {n}''''} \to {\bf n}$ corresponds to motor $i = n_{\rm pre}$ completing the first half of the stroke and motor $i = (n_{\rm pre}+n_{\rm mid})$ completing the second half of the stroke. Similarly, the release transition ${\bf n}\to{\bf{n}'}$ corresponds to the myosin $i = N$ unbinding. 

We will calculate changes to the strain values as motors transition through this FIFO {\it microstate} cycle. The strain distribution at the end of the cycle is self-consistent: it matches the strain distribution at the beginning of the cycle exactly, but with the actin position transposed and all bound motor indices incremented by one.


Beginning in {\it microstate} ${\bf {n}''}$ with $N-1$ motors bound,  we list the strain vectors for each motor in the five {\it microstate} as follows:
\begin{center}
	\resizebox{\textwidth}{!}{
	\begin{tabular}{ c  c c c c c}
		index&${\bf {n}'' }$&${\bf {n}''' }$& ${\bf {n}'''' }$& ${\bf {n}}$& ${\bf {n}'}$\\ 
		$1$&$y_{1}^{0}$ & 0 &$-\Delta x_{1}$& $-\Delta x_{1}-\Delta x_{2}$& $-\Delta x_{1}-\Delta x_{2}+\Delta x_{3}$\\ 
		\vdots&$\vdots $& $y_{1}^{0}$ & $y_{1}^{0}-\Delta x_{1}$ & $y_{1}^{0}-\Delta x_{1}-\Delta x_{2}$&  $y_{1}^{0}-\Delta x_{1}-\Delta x_{2}+\Delta x_{3}$\\  
		$ n_{\rm pre}$&$y_{n_{\rm pre}}^{0}$ & $\vdots $& $\vdots $ & $\vdots $& $\vdots $\\   
	 $n_{\rm pre+1}$+1&$y_{n_{\rm pre}+1} ^{0}+d^{\prime}$ &$y_{n_{\rm pre}}^{0}$ & $y_{n_{\rm pre}}^{0}+{\color{red}d^{\prime}}-\Delta x_{1}$  &$y_{n_{\rm pre}}^{0}+d^{\prime}-\Delta x_{1}-\Delta x_{2}$&$y_{n_{\rm pre}}^{0}+d^{\prime}-\Delta x_{1}-\Delta x_{2}+\Delta x_{3}$\\ 
		\vdots&$\vdots $&  $y_{n_{\rm pre}+1}^{0} +d^{\prime}$ & $y_{n_{\rm pre}+1}^{0} +d^{\prime}-\Delta x_{1}$&$y_{n_{\rm pre}+1}^{0} +d^{\prime}-\Delta x_{1}-\Delta x_{2}$& $y_{n_{\rm pre}+1}^{0} +d^{\prime}-\Delta x_{1}-\Delta x_{2}+\Delta x_{3}$\\ 
		$ n_{\rm pre}+n_{\rm mid}$& $	y_{n_{\rm pre}+n_{\rm mid}}^{0}+d^{\prime}$ & $\vdots $& $\vdots $& $\vdots $& $\vdots $\\ 
	$n_{\rm pre}+n_{\rm mid}+1$&$y_{n_{\rm pre}+n_{\rm mid}+1}^{0}+d$& $y_{n_{\rm pre}+n_{\rm mid}}^{0}+d^{\prime}$ &$	y_{n_{\rm pre}+n_{\rm mid}}^{0}+d^{\prime}-\Delta x_{1}$  & $	y_{n_{\rm pre}+n_{\rm mid}}^{0}+{\color{red}d}-\Delta x_{1}-\Delta x_{2}$  & $	y_{n_{\rm pre}+n_{\rm mid}}^{0}+d-\Delta x_{1}-\Delta x_{2}+\Delta x_{3}$ \\ 
	\vdots&$\vdots $ &$y_{n_{\rm pre}+n_{\rm mid}+1}^{0}+d$ &$y_{n_{\rm pre}+n_{\rm mid}+1}^{0}+d-\Delta x_{1}$&$y_{n_{\rm pre}+n_{\rm mid}+1}^{0}+d-\Delta x_{1}-\Delta x_{2}$& $y_{n_{\rm pre}+n_{\rm mid}+1}^{0}+d-\Delta x_{1}-\Delta x_{2}+\Delta x_{3}$\\ 
	$N-2$&$y_{N-2}^{0}+d $& $\vdots $& $\vdots $& $\vdots $& $\vdots $\\ 
	$N-1$&$y_{N-1}^{0}+d $& $y_{N-2}^{0}+d $& $y_{N-2}^{0}+d-\Delta x_{1} $&  $y_{N-2}^{0}+d-\Delta x_{1}-\Delta x_{2}$& $y_{N-2}^{0}+d-\Delta x_{1}-\Delta x_{2}+\Delta x_{3}$\\ 
	$N$&& $y_{N-1}^{0}+d $& $y_{N-1}^{0}+d-\Delta x_{1} $ & $y_{N-1}^{0}+d-\Delta x_{1}-\Delta x_{2} $ & 0\\ 
	\end{tabular}}
\end{center}
%\end{table}
Next, we explain the reasons for these changes in motors' strain.
The initial {\it microstate} ${\bf {n}''}$ with $N-1$ motors bound, has strain values described in Eq.~\eqref{yi_S}.
The forces exerted by these motors must balance the external load $F$ on the actin filament:
\begin{equation}
	\sum_{i=1}^{N-1} f\left(y_{i}({\bf n}'', F)\right)=F,
	\label{yi_recover}
\end{equation}
where $y_{i}({\bf n}'', F)$ is given in the first column of the table above.






 From {\it microstate} ${\bf{n}''}$ to ${\bf{n}'''}$, a new motor binds to the {\it pre-power stroke} state, and there are now $N$ motors bound. The new motor binds in a neutral position $y_{1}=0$ and all other motor indices are incremented by one, so the
 the force equilibrium equation is
 
 

\begin{equation}
	\sum_{i=1}^{N} f\left(y_{i}({\bf n}''', F)\right)=F,
	\label{yi_pre}
\end{equation}
where $y_{1}({\bf n}''', F)=0, y_{i}({\bf n}''', F)=y_{i-1}({\bf n}'', F), i=2,\ldots,N.$

From {\it microstate} ${\bf{n}'''}$ to ${\bf{n}''''}$, the $i=n_{\rm pre}$ motor changes from the {\it pre-power stroke} state to the {\it mid-power stroke} state. 
This motor strain undergoes elongation by an increment of $d^{\prime}$ nm.
 Therefore, to maintain balance with the external load $F$ on the actin filament, the strain on each motor decreases by $\Delta x_{1}({\bf n}''', F)$. The force equilibrium equation is

\begin{equation}
	\sum_{i=1}^{N} f\left(y_{i}({\bf n}'''', F)\right)=F,
	\label{yi_mid}
\end{equation}
where $y_{ n_{\rm pre}+1}({\bf n}'''', F)=y_{n_{\rm pre}+1}({\bf n}''', F)+d^{\prime}-\Delta x_{1}({\bf n}''', F),$ and $y_{i}({\bf n}'''', F)=y_{i}({\bf n}''', F)-\Delta x_{1}({\bf n}''', F)$, when $i\ne{ n_{\rm pre}+1}$.
By solving this force balance equation Eq.~\eqref{yi_mid}, we can obtain $\Delta x_{1}({\bf n}''', F)$.

Similarly, from {\it microstate} ${\bf{n}''''}$ to ${\bf n}$, the motor $i = (n_{\rm pre}+n_{\rm mid})$ transitions from {\it mid-power stroke} state into the {\it post-power stroke} state. The strain of the this motor undergoes an elongation of $d-d^{\prime}$ nm. Therefore, to maintain balance with the external load $F$ on the actin filament, the strain on each motor decreases by $\Delta x_{2}({\bf n}'''', F)$. The force equilibrium equation is

\begin{equation}
	\sum_{i=1}^{N} f\left(y_{i}({\bf n}, F)\right)=F,
	\label{yi_post}
\end{equation}
where $y_{n_{ \rm pre}+n_{\rm mid}+1}({\bf n}, F)=y_{n_{\rm pre}+n_{\rm mid}+1}({\bf n}'''', F)+(d-d^{\prime})-\Delta x_{2}({\bf n}'''', F),$ and $y_{i}({\bf n}, F)=y_{i}({\bf n}'''', F)-\Delta x_{2}({\bf n}'''', F)$, when $i\ne {n_{\rm pre}+n_{\rm mid}+1}$.
By solving this force balance equation Eq.~\eqref{yi_post}, we can obtain $\Delta x_{2}({\bf n}'''', F)$.

Finally, from {\it microstate} ${\bf n}$ to ${\bf{n}'}$, the motor $i=N$ in the {\it post-power stroke} state detaches from the actin filament.
To compensate, the $N-1$ remaining motors stretch by an amount $\Delta x_{3}({\bf n}, F) $, 
the force balance equation at this moment is
\begin{equation}
	\sum_{i=1}^{N-1} f\left(y_{i}({\bf n}', F)\right)=F,
	\label{yi_detach}
\end{equation}
where $y_{N}({\bf n}', F)=0,$ and $y_{i}({\bf n}', F)=y_{i}({\bf n}, F)+\Delta x_{3}({\bf n}, F)$, $i=1,\ldots,N-1.$
We can obtain $\Delta x_{3}({\bf n}, F)$ by Eq.~\eqref{yi_detach}. 

From {\it microstate} ${\bf {n}'}\to {\bf {n}''}$, the strain will not change as an unbound motor transitions from {\it detached} state to {\it recovered} state. The strain distribution at the end of the cycle ${\bf n}'$ return to the strain distribution at the beginning of the cycle  ${\bf n}''$. 
Let
$
	\Delta y({\bf n}'', F)=-\Delta x_{1}({\bf n}''', F)-\Delta x_{2}({\bf n}'''', F)+\Delta x_{3}({\bf n}, F)
	\label{deltay}
$, then 
\begin{equation}
	\begin{aligned}
	\Delta y({\bf n}'', F)&=y_{1}^{0}({\bf n}'', F),\\
	y_{i-1}^{0}({\bf n}'', F)+d_{i}({\bf n}'')+\Delta y({\bf n}'', F)&=y_{i}^{0}({\bf n}'', F)+d_{i}({\bf n}''), \ i=2, \ldots, N-1,
		\end{aligned}
\end{equation}
further deriving
\begin{equation}
	y_{i}^{0}({\bf n}'', F)-y_{i-1}^{0}({\bf n}'', F)=\Delta y({\bf n}'', F), \ i=2, \ldots, N-1.
\end{equation}
We find that $\left \{ y_{i}^{0}({\bf n}'', F)  \right \}_{i=1}^{N-1}$ is an arithmetic sequence with $\Delta y({\bf n}'', F)$ as the first term and $\Delta y({\bf n}'', F)$ as the common difference,
$
	y_{i}^{0}({\bf n}'', F)=i\Delta y({\bf n}'', F), \ i=1, \ldots, N-1.
$
Therefore, the strain of the motor can be expressed as
\begin{equation}
	y_{i}(\mathbf{n}'', F)=i \Delta y(\mathbf{n}'', F)+d_{i}(\mathbf{n}'').
	\label{yy_S}
\end{equation}
where $d_{i} =0,d^{\prime}$ and $d$ for a motor in the {\it pre-power stroke}, {\it mid-power stroke}, and {\it post-power stroke} states, respectively. 

The strain of $i$-th motor changes as follows:
\begin{equation}
	y_{i}= \begin{cases}i \Delta y(\mathbf{n}^{\prime \prime}, F)+d_{i}(\mathbf{n}''), & \text { in {\it microstate} } \mathbf{n}^{\prime \prime}, \\ 
	(i-1) \Delta y(\mathbf{n}^{\prime \prime}, F)+d_{i}(\mathbf{n}'''), & \text { in {\it microstate} } \mathbf{n}^{\prime \prime \prime}, \\ 
	(i-1) \Delta y(\mathbf{n}^{\prime \prime}, F)+d_{i}(\mathbf{n}'''')-\Delta x_{1}(\mathbf{n}^{\prime \prime\prime}, F), & \text { in {\it microstate} } \mathbf{n}^{\prime \prime \prime\prime}, \\
	(i-1) \Delta y(\mathbf{n}^{\prime \prime}, F)+d_{i}(\mathbf{n})-\Delta x_{1}(\mathbf{n}^{\prime \prime\prime}, F)-\Delta x_{2}(\mathbf{n}^{\prime\prime\prime \prime}, F), & \text { in {\it microstate} } \mathbf{n}, \\
	(i-1) \Delta y(\mathbf{n}^{\prime \prime}, F)+d_{i}(\mathbf{n}')-\Delta x_{1}(\mathbf{n}^{\prime\prime \prime}, F)-\Delta x_{2}(\mathbf{n}^{\prime\prime\prime \prime}, F)+\Delta x_{3}(\mathbf{n}, F), & \text { in {\it microstate} } \mathbf{n}^{\prime}, \end{cases}
\label{yii}
\end{equation}
where $\Delta y({\bf {n}''}, F)=-\Delta x_{1}({\bf {n}'''}, F)-\Delta x_{2}({\bf {n}''''}, F)+\Delta x_{3}({\bf {n}}, F)$, we can obtain $\Delta x_{1}({\bf {n}'''}, F)$, $\Delta x_{2}({\bf {n}''''}, F)$ and $\Delta x_{3}({\bf {n}}, F)$ by Eq.~\eqref{yi_mid},  Eq.~\eqref{yi_post} and Eq.~\eqref{yi_detach}. The process of changing $d_{i}$ in Eq.~\eqref{yii} is
\begin{equation}
	d_{i}= \begin{cases}d_{i}(\mathbf{n}''), & \text { in {\it microstate} } \mathbf{n}^{\prime \prime}, \\ 
		d_{i}(\mathbf{n}''')=d_{i-1}(\mathbf{n}''), & \text { in {\it microstate} } \mathbf{n}^{\prime \prime \prime}, \\ 
	d_{\rm n_{\rm pre}+1}({\bf n}'''')=d^{\prime}, d_{i}(\mathbf{n}'''')=d_{i}(\mathbf{n}'''), i\ne{\rm  n_{\rm pre}+1},& \text { in {\it microstate} } \mathbf{n}^{\prime \prime \prime\prime}, \\
	d_{\rm n_{ pre}+n_{mid}+1}({\bf n})=d, d_{i}({\bf n})=d_{i}({\bf n}''''), i\ne {\rm n_{pre}+n_{mid}+1}, & \text { in {\it microstate} } \mathbf{n}, \\
		d_{i}(\mathbf{n}')=d_{i}(\mathbf{n}''), & \text { in {\it microstate} } \mathbf{n}^{\prime}. \end{cases}
	\label{di}
\end{equation}
By calculating the last line of  Eq.~\eqref{yii}, we find that the strain in the {\it microstate}
 ${\bf{n}'}$ is $i \Delta y(\mathbf{n}^{\prime \prime}, F)+d_{i}(\mathbf{n}^{\prime \prime}) $, which is equivalent to the starting strain in {\it microstate} $\bf{{n}''}$. But, the actin filament has moved to the left a distance $-\Delta y$ and the motor indices have been incremented by 1: The motor initially indexed $i = 2$ is now indexed $i = 3$, and its strain has changed $\Delta y$, which is typically a negative quantity. This consistency across the {\it microstate}  cycle shows how our strain model of Eq.~\eqref{yy_S} captures the strain values exactly for a system with FIFO transitions.

In fact, motors transition in a first-in, first-out (FIFO) model is an overly idealized model. In this study, we employ a more generalized model that every motor in the {\it post-power stroke} state has the potential to detach. This process is stochastic, perhaps the motor with index $i=N-2$ unbinds first, rather than $i=N$. To give our model a physical interpretation, we might imagine a system in which, upon motor $i=N-2$ unbinding, the actin filament is frozen in position while the $N-1$ bound motors instantaneously rearrange to match the strain values given by our Eq.~\eqref{yy_S} for {\it microstate} ${\bf {n}''}$.



{\centering\subsection{Calculation of actin filament position}}
\begin{figure}%[tbhp]
	\centering
	\includegraphics[scale=0.27]{xn.eps}
	\caption{Schematic diagram of the change process of actin position $x_{\mathbf{n}}$. We assume that ${\bf n_{0}}=(N,0)$ is the initial position and ${{\bf {n}}=(n_{\rm pre},n_{\rm post})}$ is the final position.}
	\label{fig:1}
\end{figure}
We calculate the changes of the actin filament position $x_{\bf n}$ as it is pulled by the myosin II motors. $\Delta x_{\text{net}}$ are positive when the actin filament moves to the left. 
 In a real system, for the transition ${\bf n}\to {\bf{n}'}$, any motor in the {\it post-power storke} state can release transitions $i \in \left\lbrace n_{\text {pre}}+n_{\operatorname{mid}}+1, \ldots, N\right\rbrace $.  
 We define changes in actin position $\Delta x$ in different transitions:
\begin{equation}
	\begin{array}{ll}
		\Delta x_{\mathbf{n}^{\prime \prime}, \mathbf{n}^{\prime}}=0, & \Delta x_{\mathbf{n}^{\prime \prime \prime}, \mathbf{n}^{\prime \prime}}=0, \\
		\Delta x_{\mathbf{n}^{\prime \prime \prime \prime}, \mathbf{n}^{\prime \prime \prime}}=\Delta x_{1}\left(\mathbf{n}^{\prime \prime \prime}, F\right), & \Delta x_{\mathbf{n}, \mathbf{n}^{\prime \prime \prime \prime}}=\Delta x_{2}\left(\mathbf{n}^{\prime \prime \prime \prime}, F\right), \quad
		\Delta x_{\mathbf{n}^{\prime}, \mathbf{n}}=-\Delta x_{\text {slip}}(i, \mathbf{n}, F), 
	\end{array}
\end{equation}
where $\Delta x_{\mathbf{n}^{\prime \prime \prime \prime}, \mathbf{n}^{\prime \prime \prime}}=-\Delta x_{\mathbf{n}^{\prime \prime \prime}, \mathbf{n}^{\prime \prime \prime \prime}}$, etc.



In fact, the change in displacement $\Delta x_{\text{slip}}(i, \mathbf{n}, F)$ happens in the transitions from the {\it microstate} ${\bf n}$ to ${\bf n'}$. The slip term depends on which motor detach from actin, indexed by $i \in \left\lbrace n_{\text{pre}} + n_{\text{mid}} + 1, \ldots, N\right\rbrace $ in the {\it post-power stroke} state. To compensate, the $N-1$ remaining motors stretch by an amount $\Delta x_{\text {slip}}(i, \mathbf{n}, F) $, 
the force balance equation at this moment is
\begin{equation}
	\sum_{j=1}^{N-1} f\left(y_{j}({\bf n}', F)\right)=F,
	\label{yi_detach_i}
\end{equation}
where $y_{i}({\bf n}', F)=0,$ and $y_{j}({\bf n}', F)=y_{j}({\bf n}, F)+\Delta x_{\text {slip}}(i, \mathbf{n}, F) $, when $j\ne i.$
We can obtain $\Delta x_{\text {slip}}(i, \mathbf{n}, F) $ by Eq.~\eqref{yi_detach_i}.
As previously mentioned, when the last motor in the {\it post-power stroke} state unbinds, $\Delta x_{\text{slip}}(N, \mathbf{n}, F)$ is equal to $\Delta x_{3}({\bf n}, F)$ for $i=N$.

The average slip distance depends on the probability of any given motor releasing:
\begin{equation}
	\left\langle\Delta x_{\text {slip}}(i, \mathbf{n}, F)\right\rangle_{i}=\sum_{i=n_{\rm pre}+n_{\rm mid}+1}^{N} \Delta x_{\text {slip}}(i, \mathbf{n}, F)p(i),
\end{equation}
where $p(i)$ is the probability of  motor $i$ releasing from {\it post-power storke} state,
it depends on motors' off rates from actin, $i \in \left\lbrace n_{\rm pre }+n_{\rm mid}+1, \ldots, N\right\rbrace $,
$$p(i)=\frac{ \omega_{\text {off}}(i, \mathbf{n}, F, [\rm ATP])}{\sum_{i=n_{\rm pre}+n_{\rm mid}+1}^{N} \omega_{\rm off}(i, \mathbf{n}, F, [\rm ATP])}.$$





We define the actin position $x_{\mathbf{n}}$ of a {\it microstate}. Here $x_{\mathbf{n}}$ depends only on  $n_{\rm pre},n_{\rm mid},n_{\rm post}$, independent of $n_{\rm det},n_{\rm rec}$. 
Although we have only $n_{\rm pre}, n_{\rm post}$ axes in Fig.~\ref{fig:1}, we have the relation $n_{\rm pre}+n_{\rm mid}+n_{\rm post}=N$, $n_{\rm pre}+n_{\rm post}\le N$, for a fixed $N$. If $n_{\rm pre}$ and $n_{\rm post}$ are known, then $n_{\rm mid}$ is also known. The $(n_{\rm pre}, n_{\rm post})$ in the figure actually corresponds to $(n_{\rm pre}, n_{\rm mid}, n_{\rm post},n_{\rm det},n_{\rm rec})$, and so on in a similar fashion. There is an equivalence between coordinate transformations:
\begin{equation}
	\begin{aligned}
(n_{\rm pre}, n_{\rm post})&\Leftrightarrow (n_{\rm pre}, n_{\rm mid}, n_{\rm post},n_{\rm det},n_{\rm rec}), \\(n_{\rm pre}+1, n_{\rm post})&\Leftrightarrow(n_{\rm pre}+1, n_{\rm mid}-1, n_{\rm post},n_{\rm det},n_{\rm rec}), \\(n_{\rm pre}, n_{\rm post}-1)&\Leftrightarrow(n_{\rm pre}, n_{\rm mid}+1, n_{\rm post}-1,n_{\rm det},n_{\rm rec}).
	\end{aligned}
\end{equation}





In the Fig.~\ref{fig:1}, the distance between two adjacent points in the horizontal direction is $\Delta x_{1}({\bf n}''', F)$, and the distance between two adjacent points in the vertical direction is $\Delta x_{2}({\bf n}'''', F)$, which means that the distance between adjacent points has been obtained  in the previous section by Eq.~\eqref{yi_mid} and Eq.~\eqref{yi_post}. 

The process corresponding to the change from $(n_{\rm pre}+1, n_{\rm post})$ to $(n_{\rm pre}, n_{\rm post})$ is  the process of {\it microstate} ${\bf{n}'''}\to {\bf{n}''''}$, and 
the actin move $\Delta x_{1}({\bf n}''', F)$, where ${{\bf n}'''}=(n_{\rm pre}+1, n_{\rm mid}-1, n_{\rm post},n_{\rm det},n_{\rm rec})$. 
Similarly, the process corresponding to the change from $(n_{\rm pre}, n_{\rm post}-1)$ to $(n_{\rm pre}, n_{\rm post})$ is the process of {\it microstate} ${\bf{n}''''}\to {\bf{n}}$, and the  actin moves $\Delta x_{2}({\bf n}'''', F)$, where ${{\bf n}''''}=(n_{\rm pre}, n_{\rm mid}+1, n_{\rm post}-1,n_{\rm det},n_{\rm rec})$. 


According to Eq.~\eqref{bili}, we find that $k_{\mathbf{n}^{\prime\prime\prime\prime}, \mathbf{n}^{\prime\prime\prime}}\gg k_{\mathbf{n}^{\prime\prime\prime}, \mathbf{n}^{\prime\prime\prime\prime}}$, $k_{\mathbf{n}, \mathbf{n}^{\prime\prime\prime\prime}}\gg k_{\mathbf{n}^{\prime\prime\prime\prime}, \mathbf{n}}$, the forward rate is much larger than the reverse rate  and the forward rate is dominant.
Compared with the forward rate, the reverse rate is negligible. Therefore, we assume that the direction of motion arrow can only point upwards or to the left, as shown in the insert in Fig.~\ref{fig:1}. We search for the position of actin in a sequential order from right to left.

With some of the above explanation, we will then give the details of the method to solve for  actin position $x_{\mathbf{n}}$.
The initial position $(N,0)$ in the figure actually corresponds to ${\bf n_{0}}=(N,0,0,n_{\rm det},n_{\rm rec})$ and $x_{\bf {n_{0}}}=0$. For the convenient description, we also give the corresponding notation for the other points:
\begin{equation}
	\begin{aligned}
		(N,0)&\Leftrightarrow{\bf n_{0}}= (N,0,0,n_{\rm det},n_{\rm rec}),\\
	   (N-1, 0)&\Leftrightarrow{\bf n_{1}}=(N-1,1,0,n_{\rm det},n_{\rm rec}),\\
	   (N-1, 1)&\Leftrightarrow{\bf n_{2}}=(N-1,0,1,n_{\rm det},n_{\rm rec}),\\
	   (N-2, 0)&\Leftrightarrow{\bf n_{3}}=(N-2,2,0,n_{\rm det},n_{\rm rec}),\\
	   (N-2, 1)&\Leftrightarrow{\bf n_{4}}=(N-2,1,1,n_{\rm det},n_{\rm rec}).
	\end{aligned}
\end{equation}
The coordinate point $(N,0)$ changes to $(N-1,0)$, causing the actin filament to move a distance of $\Delta x_{1}({\bf n_{0}}, F)$, where ${{\bf n_{0}}}=(N,0,0,n_{\rm det},n_{\rm rec})$. Therefore, the actin filament's position at this stage is $x_{\mathbf{n_{1}}}=x_{\bf {n_{0}}}+\Delta x_{1}({\bf n_{0}}, F)$. Subsequently, the coordinate point changes from $(N-1,0)$ to $(N-1,1)$, resulting in a movement of the actin filament by $\Delta x_{2}({\bf n_{1}}, F)$, where ${{\bf n_{1}}}=(N-1,1,0,n_{\rm det},n_{\rm rec})$. As a result, the actin filament's position at this point is $x_{\mathbf{n_{2}}}=x_{\bf {n_{1}}}+\Delta x_{2}({\bf n_{1}}, F)$.
Similarly, the distance of change from $(N-1, 0)$ to $(N-2, 0)$ corresponds to the change in actin position, which is $\Delta x_{1}({\bf n_{1}}, F)$, where ${{\bf n_{1}}}=(N-1,1,0,n_{\rm det},n_{\rm rec})$, and the actin position is $x_{\mathbf{n_{3}}}=x_{\bf {n_{1}}}+\Delta x_{1}({\bf n}_{1}, F)$.

Next, multiple paths appear with the initial and final positions fixed.
There are two pathways from $(N, 0)$ to $(N-2, 1)$, which are {\bf path 1:} $ (N, 0)$$\to$$(N-1,0)$$\to$$(N-1,1)$$\to$$(N-2,1) $ and {\bf path 2:} $(N, 0)$$\to$$(N-1,0)$$\to$$(N-2,0)$$\to$$(N-2,1)$. 
In other words, $(N-2, 1)$ can only be changed by $(N-1,1)$ or $(N-1,0)$, which means that no matter which pathway is chosen, it must pass through one of $(N-1,1)$ or $(N-1,0)$.

If the {\bf path 1} is selected, that is, the change from $(N-1, 1)$ to $(N-2, 1)$, the distance that the actin moves is $\Delta x_{1}({\bf n}_{2}, F)$, where ${{\bf n}_{2}}=(N-1,0,1,n_{\rm det},n_{\rm rec})$.
In this {\bf path 1}, the position of  actin is $x_{\mathbf{n_{2}}}+\Delta x_{1}({\bf n}_{2}, F)$. Conversely, if the {\bf path 2} is selected, that is, the change from $(N-2, 0)$ to $(N-2, 1)$, the position of  actin is $\Delta x_{2}({\bf n}_{3}, F), {{\bf n}_{3}}=(N-2,2,0,n_{\rm det},n_{\rm rec})$.
In the {\bf path 2}, the total distance that the actin moves is $x_{\mathbf{n_{3}}}+\Delta x_{2}({\bf n}_{3}, F)$.
From the perspective of improving motor efficiency, here we choose the shorter distance between the two paths as the  position  of actin $x_{\mathbf{n_{4}}}$: 
\begin{equation}
	x_{\mathbf{n_{4}}}=\min\left\lbrace x_{\mathbf{n_{2}}}+\Delta x_{1}({\bf n}_{2}, F), x_{\mathbf{n_{3}}}+\Delta x_{2}({\bf n}_{3}, F)\right\rbrace. 
\end{equation}

And so on, according to mathematical induction,  given the final position
$x_{\mathbf{n}}=(n_{\rm pre}, n_{\rm post})$.
 Although there are many paths to change from the initial position $(N, 0)$ to the final position $(n_{\rm pre}, n_{\rm post})$, each path must pass through one of the two points $(n_{\rm pre}+1, n_{\rm post})$, $(n_{\rm pre}, n_{\rm post}-1)$. 
 If the transition from $(n_{\rm pre}+1, n_{\rm post})$ to $(n_{\rm pre}, n_{\rm post})$
 , the distance the actin moves is $\Delta x_{1}({\bf n}''', F)$, where ${{\bf n}'''}=(n_{\rm pre}+1, n_{\rm mid}-1, n_{\rm post},n_{\rm det},n_{\rm rec})$, the position of actin is $x_{\mathbf{n}'''}+\Delta x_{1}({\bf n}''', F)$.
 If the transition is made from the position $(n_{\rm pre}, n_{\rm post}-1)$ to $(n_{\rm pre}, n_{\rm post})$, the actin filament will move a distance of $\Delta x_{2}({\bf n}''''', F)$, where ${{\bf n}'''''}=(n_{\rm pre}, n_{\rm mid}+1, n_{\rm post}-1,n_{\rm det},n_{\rm rec})$. 
 The total displacement of the actin filament corresponds to the summation of its current position $x_{\mathbf{n}'''''}$, and the distance it moves, represented as $x_{\mathbf{n}} = x_{\mathbf{n}'''''} + \Delta x_{2}({\bf n}''''', F)$. To determine the precise position of the actin filament, we first calculate the distance it traverses from the initial point to the point on its right side and the point below it. And then, we incorporate the distance covered in the last step to obtain the final position.
 Finally, we then compare these two distances and choose the smaller value as the position of actin $x_{\mathbf{n}}$,
 \begin{equation}
 	x_{\mathbf{n}}=\min\left\lbrace x_{{\mathbf{n}'''}}+\Delta x_{1}({\mathbf{n}'''},F),\\ x_{\mathbf{n}''''}+\Delta x_{2}({\mathbf{n}''''},F) \right\rbrace. 
 \end{equation}
If the point lies on the hypotenuse, specifically when $n_{\rm pre} + n_{\rm post} = N$, the motor can only originate from {\bf path 2}.
  \begin{equation}
 	x_{\mathbf{n}}= x_{\mathbf{n}''''}+\Delta x_{2}({\mathbf{n}''''},F).
 \end{equation}
  By this method, we can search for the position of actin $x_{\mathbf{n}}$ at each point in Fig.~\ref{fig:1}.\\
  
{\centering\subsection{The basic free energy and transition rates}}
Here, we define the energy model that is used to calculate the microscopic transition rates and the {\it microstate} distribution $p(\mathbf{n}| N, F, {\rm [ATP]})$. We need to know the difference in basic free energy $\Delta G$ between any two {\it microstates}. 
For each transition, the exponential of the basic free energy change governs the relationship between the forward and reverse microscopic rates.
For example, 
\begin{equation}
	\frac{\omega_{\mathrm{on}}}{\omega_{\mathrm{off}}^{\prime}}=\exp \left[-\beta \Delta G_{\mathrm{bind}}\right],
	\label{G}
\end{equation}
where $\Delta G_{\rm bind}$ is the change in basic free energy upon motor binding.

There are three components to the total basic free energy. The first is {\bf the basic free energy of the individual myosin II motors}. Each motor undergoes conformational changes and hydrolyzes one ATP molecule across the cycle.  From the {\it recovered} state to the {\it pre-power stroke} state, the basic free energy of a myosin II motor changes by $g_{\rm on}$. From the {\it pre-power stroke} state and moving clockwise around the cycle, the basic free energy changes across the other transitions are labelled $g_{\rm step1}$, $g_{\rm step2}$, $g_{\rm off}$, $g_{\rm recovery}$. Due to the changes in conformational free energy sum to zero across a full cycle, the sum of these components is equal to the free energy gained from ATP hydrolysis:
\begin{equation}
	g_{\mathrm{on}}+g_{\mathrm{ps}}+g_{\mathrm{off}}+g_{\text {recovery }}=-\Delta \mu_{\mathrm{ATP}},
\end{equation}
where $\Delta \mu_{\mathrm{ATP}}\geq0$ is positive since it is a chemical potential energy of the system. The value of $\Delta \mu_{\mathrm{ATP}}$ is dependent on the concentration of ATP. For simplicity, we set $g_{\text {step1}}=g_{\text {step2} }=g_{\mathrm{ps}} / 2$. A larger absolute value of $g_{\rm ps}$ corresponds to an easier power stroke for the motor, indicating a higher likelihood of the power stroke process occurring and resulting in a greater proportion of motors transitioning from the {\it pre-power stroke} state to the {\it mid-power stroke} state and eventually to the {\it post-power stroke} state.

The second component is the {\bf strain energy of the bound motors}. The strain energy of bound myosin motors is 
\begin{equation}
	E_{\text {strain }}(\mathbf{n}, F)=\sum_{i=1}^{N}E(y_{i}({\bf{n}},F)).
	\label{EWLC_S}
\end{equation}
The third and final contribution to the basic free energy is the {\bf energy of the actin filament} as it is pulled against a constant external load, $Fx_{\bf n}$, where $x_{\bf n}$ is the change in the position of the actin filament, given in Section D.















We now calculate the local equilibrium terms $p({\bf n}|N, F, {\rm [ATP]})$ and the microscopic rates  and to solve the distribution of the system. We use this separation
of timescales in a local-equilibrium approximation that gives an analytical formula for the distribution of {\it microstates} within a given {\it mesostate}, $	p(\mathbf{n}|N, F, {\rm [ATP]})$. 



First, the local equilibrium approximation solves for the conditional distribution of all {\it microstates} within a given {\it mesostate} $N$. This {\it microstate} conditional distribution contains a product of two terms that correspond to relaxation within the bound motors and within the unbound motors \cite{wagoner2021evolution}:

\begin{equation}
	\begin{aligned}
		p(\mathbf{n}|N, F, {\rm [ATP]})=& \frac{N !}{n_{\rm pre } ! n_{\rm mid} ! n_{\rm post} !} \frac{\exp \left[-\beta G_{B}(\mathbf{n}, F, {\rm [ATP]})\right]}{Q_{B}}\times \frac{\left(N_{T}-N\right) !}{n_{\text {det}} ! n_{\text {rec}} !} \frac{\exp \left[-\beta G_{U}(\mathbf{n},{\rm [ATP]})\right]}{Q_{U}},
		\label{pn}
	\end{aligned}
\end{equation}
where $\beta=1/k_{B}T$, and
\begin{equation}
	\begin{aligned}
		G_{B}(\mathbf{n}, F, {\rm [ATP]})=&\left[\frac{n_{\rm mid}}{2}+n_{\rm post}\right]  g_{\mathrm{ps}}+E_{\text {strain }}\left(n_{\rm pre}, n_{\rm mid}, n_{\rm post}, F\right)-E_{\text {strain }}\left(N, 0, 0, F\right)+Fx_{\bf{n}},
	\end{aligned}
\label{GB}
\end{equation}
is the energy of the bound motors relative to the {\it microstate} in which all bound motors are in the {\it pre-power stroke} state. Similarly,
\begin{equation}
	G_{U}(\mathbf{n},{\rm [ATP]})=n_{\text {rec}} g_{\text {recovery}},
\end{equation}
is the energy of the unbound motors relative to the {\it microstate} in which all unbound motors are in the {\it detached} state. $Q_{B}$ and $Q_{U}$ are the normalization factors for these two terms. The following two equations are satisfied:
\begin{equation}
	\begin{aligned}
&\sum_{n_{\rm pre}+n_{\rm mid}+n_{\rm post}=N}^{} \frac{N !}{n_{\rm pre} ! n_{\rm mid} ! n_{\rm post} !} \frac{\exp \left[-\beta G_{B}(\mathbf{n}, F, {\rm [ATP]})\right]}{Q_{B}} =1,\\
&\sum_{n_{\text {det}}+n_{\text {rec}}=N_{T}-N}^{}\frac{\left(N_{T}-N\right) !}{n_{\text {det}} ! n_{\text {rec}} !} \frac{\exp \left[-\beta G_{U}(\mathbf{n}, {\rm [ATP]})\right]}{Q_{U}}=1.
	\end{aligned}
\end{equation}









We introduce the energy change upon motor binding to the {\it pre-power stroke} state and upon motor release from the {\it post-power stroke} state:
\begin{equation}
	\begin{aligned}
		&\Delta G_{\rm bind}\left(\mathbf{n}^{\prime \prime}, F, {\rm [ATP]}\right)= g_{\mathrm{on}}+E_{\text {strain }}\left(\mathbf{n}^{\prime \prime \prime}, F\right)-E_{\text {strain }}\left(\mathbf{n}^{\prime \prime}, F\right), \\
		&\Delta G_{\rm rel}(i, \mathbf{n}, F, {\rm [ATP]})= g_{\text {off}}+E_{\text {strain }}\left(\mathbf{n}^{\prime}, F\right)-E_{\text {strain }}(\mathbf{n}, F)-F\Delta x_{\text {slip}}(i, \mathbf{n}, F).
	\end{aligned}
\end{equation}
 

For the transition between $\mathbf{n}^{\prime\prime}$ and $\mathbf{n}^{\prime\prime\prime}$, the attachment rate of myosin with actin filament, $\omega_{\mathrm{on}}\left(\mathbf{n}^{\prime \prime}, {\rm [ATP]}\right)$ and $\omega_{\mathrm{off}}^{\prime}\left(\mathbf{n}^{\prime \prime \prime}, F\right)$, is govern by the Boltzmann equation as follows:
\begin{equation}
	\begin{aligned}
		&\omega_{\mathrm{on}}\left(\mathbf{n}^{\prime \prime}, {\rm [ATP]}\right) =\omega_{\mathrm{on}}^{0} e^{-\beta g_{\mathrm{on}}}, \\
		&\omega_{\mathrm{off}}^{\prime}\left(\mathbf{n}^{\prime \prime \prime}, F\right) =\omega_{\mathrm{on}}^{0} e^{\beta\left(\Delta G_{\mathrm{bind}}\left(\mathbf{n}^{\prime \prime}, F\right)-g_{\mathrm{on}}\right)},
	\end{aligned}
\end{equation}
where $\omega_{\rm on}^{0}$ is the intrinsic rate.

The {\it post-power stroke} rates depend on not just the {\it microstate}, but also the index of the motor $i$,
\begin{equation}
	\begin{aligned}
		\omega_{\mathrm{off}}(i, \mathbf{n}, F, {[\rm ATP]})=& \begin{cases}{[\rm ATP]}\omega_{\mathrm{off}}^{0} e^{-\beta \delta f(y_{i}(\mathbf{n}, F))}, & i \in\left\lbrace n_{\rm pre}+n_{\rm mid}+1, \ldots, N\right\rbrace , \\
			0 ,& \text {otherwise},\end{cases} \\
		\omega_{\mathrm{on}}^{\prime}\left(i, \mathbf{n}^{\prime}, F, {[\rm ATP]} \right)=& \begin{cases}\omega_{\mathrm{off}}^{0} e^{-\beta \delta f(y_{i}(\mathbf{n}, F))} e^{\beta \Delta G_{\mathrm{rel}}(i,\mathbf{n}, F)}, & i \in\left\lbrace n_{\rm pre}+n_{\rm mid}+1, \ldots, N\right\rbrace , \\
			0, & \text {otherwise},\end{cases}
	\end{aligned}
\label{omega_off_on}
\end{equation}
where $\omega_{\rm off}^{0}$ is the off-rate at zero force. The rate $\omega_{\mathrm{off}}$ and $\omega_{\mathrm{on}}^{\prime}$ are related to the external load $F$ and ATP concentration [ATP].
For the rate $\omega_{\rm on}^{\prime}$: to ensure reversibility, the system begins in ${\bf {n}^{\prime}}$ and transitions to $\bf n$ when a motor in the {\it detached} state spontaneously stretches the appropriate amount to bind in a {\it post-power stroke}
position with force $f(y_{i}({\bf n},F))$, where $i$ gives the index over {\it post-power stroke} motors in the {\it microstate} $\bf n$.




Next, we give the transition rate between {\it pre-power stroke} state and {\it mid-power stroke} state  and the transition rate between {\it mid-power stroke} state  and {\it post-power stroke} state:

\begin{equation}
	\begin{array}{ll}
		\omega_{\rm pre}\left(\mathbf{n}^{\prime \prime\prime}, F, {\rm [ATP]}\right)=\omega_{\rm pre}^{0}e^{-\beta\left[ \Delta G_{1}\left(\mathbf{n}^{\prime \prime\prime}, F\right)+g_{\mathrm{ps}}/2\right] }, & \omega^{\prime}_{\rm pre}\left(\mathbf{n}^{\prime \prime\prime\prime}, {\rm [ATP]}\right)= \omega_{\rm pre}^{0}e^{-\beta g_{\mathrm{ps}}/2}, \\
		\omega_{\rm mid}\left(\mathbf{n}^{\prime \prime\prime\prime}, F, {\rm [ATP]}\right)= \omega_{\rm mid}^{0}e^{-\beta\left[ \Delta G_{2}\left(\mathbf{n}^{\prime\prime \prime\prime}, F\right)+g_{\mathrm{ps}}/2\right] }, &  \omega^{\prime}_{\rm mid}\left(\mathbf{n}, {\rm [ATP]}\right)= \omega_{\rm mid}^{0}e^{-\beta g_{\mathrm{ps}}/2},
	\end{array}
\end{equation}
where $\omega_{\rm pre}^{0}$ and $\omega_{\rm mid}^{0}$ are load free transition rates.

We need the energy change from  {\it pre-power stroke} to  {\it mid-power stroke} and from  {\it mid-power stroke} to {\it post-power stroke}:   
\begin{equation}
	\begin{aligned}
		\Delta G_{1}\left(\mathbf{n}^{\prime \prime\prime}, F, {\rm [ATP]}\right)=& \frac{g_{\mathrm{ps}}}{2}+E_{\text {strain }}\left(\mathbf{n}^{\prime\prime \prime \prime}, F\right)-E_{\text {strain }}\left(\mathbf{n}^{\prime\prime \prime}, F\right)+F\Delta x_{1}(\mathbf{n}^{\prime\prime \prime}, F), \\
		\Delta G_{2}(\mathbf{n}^{\prime \prime\prime\prime}, F, {\rm [ATP]})=& \frac{g_{\mathrm{ps}}}{2}+E_{\text {strain }}\left(\mathbf{n}, F\right)-E_{\text {strain }}(\mathbf{n}^{\prime \prime\prime\prime}, F)+F\Delta x_{2}(\mathbf{n}^{\prime\prime\prime\prime}, F),
	\end{aligned}
\end{equation}
where we are once again using notation from the {\it microstate} cycle, Eq.~\eqref{n5_S}, for convenience.


Similar to the relationship described in Eq.~\eqref{G}, the exponential of the basic free energy change constrains the ratio of forward and reverse microscopic rates. The transition from the {\it pre-power stroke} state to the {\it mid-power stroke} state and from the {\it mid-power stroke} state to the {\it post-power stroke} state also follows a rate relationship, which can be derived as follows:
\begin{equation}
	\begin{aligned}
		\frac{\omega_{\rm pre}}{\omega^{\prime}_{\rm pre}}=e^{-\beta \Delta G_{1}\left(\mathbf{n}^{\prime \prime\prime}, F\right)},\quad\frac{\omega_{\rm mid}}{\omega^{\prime}_{\rm mid}}=e^{-\beta \Delta G_{2}\left(\mathbf{n}^{\prime\prime \prime\prime}, F\right)},
	\end{aligned}
\label{bili}
\end{equation}
where $\omega_{\rm pre}$ and $\omega_{\rm mid}$ are microscopic transition rates from  {\it pre-power stroke} state to  {\it mid-power stroke} state and from  {\it mid-power stroke} state to {\it post-power stroke} state of an individual motor. $\omega^{\prime}_{\rm pre}$ is the reverse rate corresponding to $\omega_{\rm pre}$, 
and $\omega^{\prime}_{\rm mid}$
is the reverse rate corresponding to $\omega_{\rm mid}$.






According to the FIFO manner we proposed earlier, only one motor has undergone a state change from {\it pre-power stroke} state to {\it mid-power stroke} state, and only one motor changes from {\it mid-power stroke} state to {\it post-power stroke} state, so microscopic rates of transfer between $\mathbf{n}^{\prime\prime\prime}$ and $\mathbf{n}^{\prime\prime\prime\prime}$ and between $\mathbf{n}^{\prime\prime\prime\prime}$ and $\mathbf{n}$ are
\begin{equation}
	\begin{array}{ll}
		k_{\mathbf{n}^{\prime\prime\prime\prime}, \mathbf{n}^{\prime\prime\prime}}=	\omega_{\rm pre}\left(\mathbf{n}^{\prime \prime \prime}, F, {\rm [ATP]}\right), & k_{\mathbf{n}^{\prime\prime \prime}, \mathbf{n}^{\prime\prime \prime \prime}}= \omega^{\prime}_{\rm pre}\left(\mathbf{n}^{\prime \prime \prime\prime}, {\rm [ATP]}\right), \\
		k_{\mathbf{n}, \mathbf{n}^{\prime\prime\prime\prime}}=\omega_{\rm mid}\left(\mathbf{n}^{\prime \prime \prime\prime}, F, {\rm [ATP]}\right), & k_{\mathbf{n}^{\prime \prime\prime \prime}, \mathbf{n}}=\omega^{\prime}_{\rm mid}\left(\mathbf{n}, {\rm [ATP]}\right).\\
	\end{array}
\end{equation}\\





{\centering\subsection {Calculating the steady state mesostate distribution}}
The mesoscopic rates $k_{N+1, N},k_{N-1, N}^{\prime},k_{N-1, N},k_{N+1, N}^{\prime}$ of linear network can be written as

\begin{equation}
	\begin{aligned}
		k_{N+1, N}&=\sum_{\mathbf{n}} p\left(\mathbf{n}^{\prime \prime}|N, F, {\rm [ATP]}\right)n_{\text {rec}} \omega_{\mathrm{on}}\left(\mathbf{n}^{\prime \prime}, {\rm [ATP]}\right),\\
		k_{N-1, N}^{\prime} &=\sum_{\mathbf{n}} p\left( \mathbf{n}^{\prime \prime\prime}|N, F, {\rm [ATP]}\right) n_{\text {pre}} \omega_{\mathrm{off}}^{\prime}\left(\mathbf{n}^{\prime \prime \prime}, F\right), \\
		k_{N-1, N} &=\sum_{\mathbf{n}} p\left( \mathbf{n}|N, F, {\rm [ATP]}\right) \sum_{i=n_{\rm pre}+n_{\rm mid}+1}^{N} \omega_{\mathrm{off}}(i, \mathbf{n}, F, [\rm ATP]) , \\
		k_{N+1, N}^{\prime} &=\sum_{\mathbf{n}} p\left( \mathbf{n}^{\prime}|N, F, {\rm [ATP]}\right) n_{\det}\sum_{k=n_{\rm pre}+n_{\rm mid}+1}^{N+1} \omega_{\text {on}}^{\prime}\left(k, \mathbf{n}^{\prime}, F, [\rm ATP]\right) p^{\prime}(k).
	\end{aligned}
\end{equation}

\begin{figure}%[tbhp]
	\centering
	\includegraphics[scale=0.25]{markovchain.eps}
	\caption{The linear network of {\it mesostates} of the myosin II motors system. We project a system of $N_{\rm T}$ motors onto this network of states, in which $N$ is the number of bound motors.}
	\label{fig:2}
\end{figure}
We solve for the steady state distribution of an actomyosin filament at fixed
external load $F$ and [ATP], which is experimentally controllable. We now calculate the {\it mesostate} distribution $p(N|F, [\rm ATP])$. This is a linear Markov chain of {\it mesostates} rates. 
The forward rate from $N$ to $N+1$ is $g\left( N\right) $, the reverse rate from $N+1$ to $N$ is $r\left( N+1\right)$, $N\in \left\lbrace 1,2,\ldots,N_{T}-1\right\rbrace $,
\begin{equation}
\begin{aligned}
		&g(N-1)=k_{N, N-1}+k_{N, N-1}^{\prime}, && g(N)=k_{N+1, N}+k_{N+1, N}^{\prime}, \\
		&r(N)=k_{N-1, N}+k_{N-1,N}^{\prime}, && r(N+1)=k_{N, N+1}+k_{N, N+1}^{\prime}.
\end{aligned}
\end{equation}
The probability $p_{N}(F, [{\rm ATP}], t)$ that $N$ motors are bound at time $t$ obeys the one-step master equation
\begin{equation}
	\frac{dp_{N}}{d t} (F, [{\rm ATP}], t)=g(N-1)p_{N-1}(F, [{\rm ATP}], t)+r(N+1)p_{N+1}(F, [{\rm ATP}], t) -[r(N)+g(N)]p_{N}(F, [{\rm ATP}], t).
\end{equation}
Let $p(N|F,{\rm [ATP]})=\underset{t\to\infty}{\lim}p_{N}(F, {\rm [ATP]}, t)$ be the steady-state {\it mesostate}  probability of $N$ motors, we have for $2\le N \le N_{T}$,
\begin{equation}
g(N-1) p(N-1|F, {\rm [ATP]})+r(N+1)p(N+1|F, {\rm [ATP]})-[r(N)+g(N)] p(N|F, {\rm [ATP]})=0.
\label{steady}
\end{equation}
In particular, the first equation is
\begin{equation}
r(2) p(2|F, {\rm [ATP]})-g(1) p(1|F, {\rm [ATP]})=0.
\end{equation}
From Eq.~\eqref{steady}, we derive that
\begin{equation}
	\begin{aligned}
		r(N+1) p(N+1|F, {\rm [ATP]})-g(N) p(N|F, {\rm [ATP]}) &=r(N) p(N|F, {\rm [ATP]})-g(N-1) p(N-1|F, {\rm [ATP]}) \\
		&=\cdots \\
		&=r(2) p(2|F, {\rm [ATP]})-g(1) p(1|F, {\rm [ATP]}) \\
		&=0.
	\end{aligned}
\end{equation}
We obtain the recursive formula for the probability
\begin{equation}
		\frac{p(N|F, {\rm [ATP]})}{p(N-1|F, {\rm [ATP]})}=\frac{g(N-1)}{r(N)} .
\end{equation}
Further deducing that
\begin{equation}
		\frac{p(N|F, {\rm [ATP]})}{p(1|F, {\rm [ATP]})}=\frac{p(N|F, {\rm [ATP]})}{p(N-1|F, {\rm [ATP]})} \cdot \frac{p(N-1|F, {\rm [ATP]})}{p(N-2|F, {\rm [ATP]})} \cdots \frac{p(2|F, {\rm [ATP]})}{p(1|F, {\rm [ATP]})}=\frac{g(N-1)}{r(N)} \cdot \frac{g(N-2)}{r(N-1)} \cdots \frac{g(1)}{r(2)}.
\end{equation}
The simplification gives
\begin{equation}
	p(N|F, {\rm [ATP]})=\prod_{j=1}^{N-1} \frac{g(j)}{r(j+1)} p(1|F, {\rm [ATP]}),
	\label{pNN}
\end{equation}
and then $p(1|F, {\rm [ATP]})$ can be obtained from the normalization condition
\begin{equation}
	\sum_{k=1}^{N_{T}} p(k|F, {\rm [ATP]})=1,
\end{equation}
which implies
\begin{equation}
	\left(1+\sum_{k=2}^{N_{T}} \prod_{j=1}^{k-1} \frac{g(j)}{r(j+1)}\right) p(1|F, {\rm [ATP]})=1.
\end{equation}
Therefore, we obtain that
\begin{equation}
	p(1|F, {\rm [ATP]})=\frac{1}{1+\sum_{k=2}^{N_{T}} \prod_{j=1}^{k-1} \frac{g(j)}{r(j+1)}}.
	\label{p0}
\end{equation}
Combining Eq.~\eqref{pNN} and Eq.~\eqref{p0}, the steady-state {\it mesostate} probability  is
\begin{equation}
	p(N|F, {\rm [ATP]})=\frac{\prod_{j=1}^{N-1} \frac{g(j)}{r(j+1)}}{1+\sum_{k=2}^{N_{T}} \prod_{j=1}^{k-1} \frac{g(j)}{r(j+1)}}=\frac{\prod_{j=1}^{N-1} \frac{k_{j+1, j}+k_{j+1, j}^{\prime}}{k_{j, j+1}+k_{j, j+1}^{\prime}}}{1+\sum_{k=2}^{N_{T}} \prod_{j=1}^{k-1} \frac{k_{j+1, j}+k_{j+1, j}^{\prime}}{k_{j, j+1}+k_{j, j+1}^{\prime}}}.
	\label{pN}
\end{equation}
Therefore, we solve for $p(N|F, {\rm [ATP]})$ using the known analytical solution for linear Markov chains. Then, combining Eq.~\eqref{pn} and Eq.~\eqref{pN}, we can obtain 
\begin{equation}
p ({\bf n}|F, {\rm [ATP]}) = p ({\bf n},N|F, {\rm [ATP]}) = p({\bf n}|N, F, {\rm [ATP]})p(N|F, {\rm [ATP]}).
\end{equation}
The three-dimensional plot in Fig.~\ref{3Dprobability} displays the steady-state {\it mesostate} probability $p(N|F)$ and the conditional distribution $p({\bf n}|N,F)$, providing a more intuitive visualization.
\begin{figure}[htbp]
	\includegraphics[scale=0.45]{3DP.eps}\\
	\caption{\label{3Dprobability} 3D diagram of the steady-state {\it mesostate} probability $p(N|F)$ and the {\it microstate} conditional distribution $p({\bf n}|N,F)$. In this figure, [ATP]=2000 $\mu$m. \textbf{(a)} The steady-state {\it mesostate} probability $p(N|F)$ is plotted against the number of bound motors $N$ and the load $F$.
	\textbf{(b-d)} The {\it microstate} conditional distribution $p({\bf n}|N, F)$ versus the number of motors in {\it pre-power stroke} state $n_{\rm pre}$ and {\it post-power stroke} state $n_{\rm post}$.}
\end{figure}	\\

{\centering\subsection{The expression of average number of bound motors, velocity  and sliding distance}}

We can define the average number of bound motors as
\begin{equation}
	\langle N(F, {[\rm ATP]})\rangle=\sum_{N=1}^{N_{T}} Np(N|F, {\rm [ATP]}).
\end{equation}

Next, we could define the velocity of the actin filament as
 \begin{equation}
	\begin{aligned}
		V(F, [\rm ATP])&=\sum_{\mathbf{n}}\Big[ p\left(\mathbf{n}^{\prime\prime\prime}\mid F, {\rm [ATP]}\right)k_{\mathbf{n}^{\prime\prime\prime\prime}, \mathbf{n}^{\prime\prime\prime}}\Delta x_{1}\left(\mathbf{n}^{\prime \prime \prime}, F\right)
		+p\left(\mathbf{n}^{\prime\prime\prime\prime}\mid F, {\rm [ATP]}\right)k_{\mathbf{n}^{\prime\prime \prime}, \mathbf{n}^{\prime\prime \prime \prime}}\left( -\Delta x_{1}\left(\mathbf{n}^{\prime \prime \prime}, F\right)\right)\Big]  \\
		&\quad+\sum_{\mathbf{n}}\Big[\left(p\left(\mathbf{n}^{\prime\prime\prime\prime}\mid F, {\rm [ATP]}\right)k_{\mathbf{n}, \mathbf{n}^{\prime\prime\prime\prime}}\Delta x_{2}\left(\mathbf{n}^{\prime \prime\prime \prime}, F\right)+p\left(\mathbf{n}\mid F, {\rm [ATP]}\right)k_{\mathbf{n}^{\prime\prime\prime\prime}, \mathbf{n}}\left( -\Delta x_{2}\left(\mathbf{n}^{\prime\prime \prime \prime}, F\right)\right)\right)\Big]  \\
		&\quad+\sum_{\mathbf{n}}\Big[
		p\left(\mathbf{n}\mid F, {\rm [ATP]}\right)\sum_{i=n_\text{pre}+n_\text{mid}+1}^{N} \omega_{\mathrm{off}}\left( i, \mathbf{n}, F, {\rm [ATP]})\left( -\Delta x_{\text {slip}}(i, \mathbf{n}, F\right) \right) p(i)\Big]\\&	 \quad+\sum_{\mathbf{n}}\Big[ p\left(\mathbf{n}^{\prime}\mid F, {\rm [ATP]}\right)\sum_{k=n_\text{pre}+n_\text{mid}+1}^{N} \omega_{\text {on}}^{\prime}\left(k, \mathbf{n}^{\prime}, F, {\rm [ATP]}\right) \Delta x_{\text {slip}}(k, \mathbf{n}, F)p^{\prime}(k)\Big],
		\label{SV}
	\end{aligned}
\end{equation}
where $p(i)$ is the probability of  motor $i$ releasing from {\it post-power storke} state,
which depends on the off rates of motors from actin and can be represented as
$$p(i)=\frac{ \omega_{\rm off }(i, \mathbf{n}, F, [\rm ATP])}{\sum_{i=n_\text{pre}+n_\text{mid}+1}^{N} \omega_{\rm off }(i, \mathbf{n}, F, [\rm ATP])},$$
and $p^\prime(k)$ is the probability of a detached motor spontaneously stretching the appropriate amount to bind in a {\it post-power stroke} position $k$ with force $f(y_{k})$, where $k$ gives the index over {\it post-power stroke} motors in the final {\it microstate} ${\bf n}$. The value of $p^\prime(k)$ depends on the attachment rates of motors,
$$p^{\prime}(k)=\frac{ \omega_{\rm on}^{\prime}(k, \mathbf{n}^{\prime}, F, [\rm ATP])}{\sum_{k=n_\text{pre}+n_\text{mid}+1}^{N} \omega_{\rm on}^{\prime}(k, \mathbf{n}^{\prime}, F, [\rm ATP])}.$$
 

 
 \begin{figure}%[tbhp]
 	\centering
 	\includegraphics[scale=0.55]{S1.eps}
 	\caption{Plot of the three components of $\left\langle\Delta x_{\rm net}\right\rangle$  and $L$. \textbf{(a)} Division of $\left\langle\Delta x_{\rm net}\right\rangle$  into three distinct components, and their corresponding average distances for the three processes in the cycle: ${\bf {n}'''}\to {\bf {n}''''}$ ($\left\langle\Delta x_{1}\right\rangle$), ${\bf {n}''''\to {n}}$ ($\left\langle\Delta x_{2}\right\rangle$), and ${\bf n}\to {\bf{n}'}$ ($\left\langle\Delta x_{\rm slip}\right\rangle$). \textbf{(b)} Similarly, we divided $L$ into three distinct components to analyze the contribution of each sub-process. }
 	\label{threepartsXnet}
 \end{figure}
 
 In addition, $\Delta x_{\text{net}}$ denotes the net distance of each step taken by the actin filament in one  motor cycle. The average net distance can be defined as
  \begin{equation}
 	\begin{aligned}
 		\left\langle\Delta x_{\rm net}(F, {[\rm ATP]})\right\rangle&=\sum_{\mathbf{n}}\Big[  \frac {p\left(\mathbf{n}^{\prime\prime\prime}\mid F, {\rm [ATP]}\right)}{\sum_{\mathbf{n}}p\left(\mathbf{n}^{\prime\prime\prime}\mid F, {\rm [ATP]}\right)}\Delta x_{1}\left(\mathbf{n}^{\prime \prime \prime}, F\right) +\frac{ p\left(\mathbf{n}^{\prime\prime\prime\prime}\mid F, {\rm [ATP]}\right)}{\sum_{\mathbf{n}}p\left(\mathbf{n}^{\prime\prime\prime\prime}\mid F, {\rm [ATP]}\right)}\Delta x_{2}\left(\mathbf{n}^{\prime\prime \prime \prime}, F\right)\Big] \\
 		&\quad+\sum_{\mathbf{n}}\Big[ \frac{p\left(\mathbf{n}\mid F, {\rm [ATP]}\right)}{\sum_{\mathbf{n}}p\left(\mathbf{n}\mid F, {\rm [ATP]}\right)}\sum_{i=n_\text{pre}+n_\text{mid}+1}^{N}\left( -\Delta x_{\text {slip}}(i, \mathbf{n}, F)\right) p(i)\Big] .
 	\end{aligned}
 \label{Xnet}
 \end{equation}
To further analyze $\left\langle\Delta x_{\rm net}(F, {[\rm ATP]})\right\rangle$ based on the cycle changes across ${\bf {n}'''}\to {\bf {n}''''}$ (step 1), ${\bf {n}''''\to {n}}$ (step 2), and ${\bf n}\to {\bf{n}'}$ (slip), we can divide it into three distinct components and calculate the average distance for each of these processes separately.
  \begin{equation}
 	\begin{aligned}
 		&\left\langle\Delta x_{ 1}(F, {[\rm ATP]})\right\rangle=\sum_{\mathbf{n}}\Big[ \frac {p\left(\mathbf{n}^{\prime\prime\prime}\mid F, {\rm [ATP]}\right)}{\sum_{\mathbf{n}}p\left(\mathbf{n}^{\prime\prime\prime}\mid F, {\rm [ATP]}\right)}\Delta x_{1}\left(\mathbf{n}^{\prime \prime \prime}, F\right) \Big],\\
 		&\left\langle\Delta x_{2}(F, {[\rm ATP]})\right\rangle=\sum_{\mathbf{n}}\Big[ \frac{ p\left(\mathbf{n}^{\prime\prime\prime\prime}\mid F, {\rm [ATP]}\right)}{\sum_{\mathbf{n}}p\left(\mathbf{n}^{\prime\prime\prime\prime}\mid F, {\rm [ATP]}\right)}\Delta x_{2}\left(\mathbf{n}^{\prime\prime \prime \prime}, F\right)\Big],\\
 		&\left\langle\Delta x_{\rm slip}(F, {[\rm ATP]})\right\rangle=\sum_{\mathbf{n}}\Big[ \frac{p\left(\mathbf{n}\mid F, {\rm [ATP]}\right)}{\sum_{\mathbf{n}}p\left(\mathbf{n}\mid F, {\rm [ATP]}\right)}\sum_{i=n_\text{pre}+n_\text{mid}+1}^{N}\left( -\Delta x_{\text {slip}}(i, \mathbf{n}, F)\right) p(i)\Big].
 	\end{aligned}
 	\label{threepartsofXnet}
 \end{equation}
 Next, we plot each of these three components of $\left\langle\Delta x_{\rm net}(F, {[\rm ATP]})\right\rangle$ separately, as illustrated in the Fig.~\ref{threepartsXnet}(a), enabling us to clearly visualize the contributions of each sub-process to $\left\langle\Delta x_{\rm net}(F, {[\rm ATP]})\right\rangle$ and how they vary in response to the load $F$.

We also obtain the  sliding distance $L$. The work per cycle of motor  is $F\Delta x_{\rm net}$, the sliding distance $L$ is calculated as work divided by the force per motor $F/N$. This motivates the relation $L=\left\langle N\Delta x_{\rm net}\right\rangle$.
In a similar way, we can divide $L$ into three distinct components and plot their corresponding graphs in Fig.~\ref{threepartsXnet}(b).
  \begin{equation}
	 L_{1}=\left\langle N\Delta x_{1}(F, {[\rm ATP]})\right\rangle,\quad
	 L_{2}=\left\langle N\Delta x_{2}(F, {[\rm ATP]})\right\rangle, \quad
	 L_{\rm slip}=\left\langle N\Delta x_{\rm slip}(F, {[\rm ATP]})\right\rangle.
	\label{threepartsofL}
\end{equation}
The thermodynamic efficiency of the myosin half-filament is
\begin{equation}
	\eta(F, {[\rm ATP]})=\frac{F\left\langle\Delta x_{\mathrm{net}}(F, {[\rm ATP]})\right\rangle}{\Delta \mu_{\mathrm{ATP}}},
\end{equation}
where $F\left\langle \Delta x_{\text {net}}\right\rangle$ is the average work done on the actin filament for one full motor cycle and $\Delta \mu_{\rm ATP}\ge 0$ is the free energy gained from ATP hydrolysis across a full cycle.


The free energy change $\Delta \mu_{\rm ATP}$ resulting from ATP hydrolysis is dependent on the concentration of ATP. The equation governing ATP hydrolysis is given by
\begin{equation}
	\Delta \mu_{\rm ATP}=\Delta \mu_{\rm ATP}^{0}+k_{B}T\ln{[\rm ATP]},
	\label{ATPP}
\end{equation}
where $\Delta \mu_{\rm ATP}^{0}$ is the ratio of concentrations from the standard value. Under standard conditions, where $[\rm ATP]=2000\mu M$, the free energy change of ATP hydrolysis is 25 $k_{ B}T$. We take $k_{B}T$ as the unit of energy. Thus, the above equation becomes
\begin{equation}
	25 \ k_{ B}T=\Delta \mu_{\rm ATP}^{0}+k_{B}T\ln{2000}.
		\label{ATP2000}
\end{equation}
By combining Eqs. \eqref{ATPP} and \eqref{ATP2000}, we obtain the following equation
\begin{equation}
	\Delta \mu_{\rm ATP}=\left( 25 +\ln{\frac{[\rm ATP]}{2000}}\right) \ k_{ B}T .
	\label{ATPG_S}
\end{equation}
As seen in Eq. \eqref{ATPG_S},  $\Delta\mu_{\rm ATP}$ increases with an increase in ATP concentration.

\begin{figure}[htbp]
	\includegraphics[scale=0.5]{ATP.eps}
	\caption{\label{ATP} 
		Biophysical properties of myosin II motors as functions of ATP concentration under different loads. \textbf{(a)} Variation of velocity $V$. \textbf{(b)} Variation of average number of bound motors $\left\langle N\right\rangle$. \textbf{(c)} Variation of force per motor $F/\left\langle N\right\rangle $. \textbf{(d)} Variation of sliding distance $L$. \textbf{(e)} Variation of net distance $\left\langle \Delta x_{\text{net}}\right\rangle $. \textbf{(f)} Variation of output work $W$. When the ATP concentration changes, we maintain the proportion of energy allocated to each component during ATP hydrolysis constant, that is, we keep $g_{\rm on}:g_{\rm off}:g_{\rm ps}:g_{\rm recover}$ constant.}
\end{figure}

\begin{figure}[htbp]
	\includegraphics[scale=0.5]{Figure1_F_ATP.eps}\\
	\caption{\label{Figure1_F_ATP}This figure serves as an extension of the findings presented in Fig.~3. Here, we illustrate the variations of biophysical properties with respect to $F$ at different ATP concentrations. In contrast, Fig.~3 specifically illustrates the fitting results and experimental data at the physiological ATP concentration of [ATP]=2000 $\mu$M.}
\end{figure}


{\centering\subsection{The mean run time/length, mean existence probability and half-life period of myosin filament}}





Let $F_N(t)$ be the probability density function representing the first-time separation of myosin filament from the actin filament (i.e., reaching 0 motor binding) at time $t$, with $N$ motors initially binding at $t=0$, $N\in \left\lbrace 1,2,\ldots,N_{T}\right\rbrace $.
This enables us to investigate the behavior of myosin II motors on the actin filament. It can be shown that $F_N(t)$ satisfy the following backward master equations, 

\begin{equation}
	\begin{aligned}
		\frac{d F_{1}(t)}{d t}=&r(1) F_{0}(t)+g(1) F_{2}(t)-[r(1)+g(1)] F_{1}(t),\\
		&\cdots \\
		\frac{d F_{N(t)}}{d t}=&r(N) F_{N-1}(t)+g(N) F_{N+1}(t)-[r(N)+g(N)] F_{N}(t), \\
		&\cdots \\
		\frac{d F_{N_{T}}}{d t}=&r\left(N_{T}\right) F_{N_{T}-1}(t)-r\left(N_{T}\right) F_{N_{T}}(t),
		\label{T}
	\end{aligned}
\end{equation}
where $F_{0}(t)=\delta(t)$ in the first equation means that if myosin filament detach from actin at $t=0$, the first-passage process is accomplished immediately, and $r(1)$ is the detachment rate of myosin filament from actin.

The run time of myosin filament initiated from $N$ motors bound is
\begin{equation}
	T_{N}=\int_{0}^{+\infty} t F_{N}(t) d t.
\end{equation}
By multiplying $t$ and integrating both sides of Eq.~\eqref{T}, the left side of the equation is
\begin{equation}
	\int_{0}^{+\infty} t \frac{d F_{N}(t)}{d t} d t=t F_{N}(t)\Big|_{0} ^{+\infty}-\int_{0}^{+\infty} F_{N}(t) d t=-1,
\end{equation}
and the first term on the right-hand side of the first equation is
\begin{equation}
	\int_{0}^{+\infty} t \delta(t) d t=0.
\end{equation}
Then, from Eq.~\eqref{T}, we have

\begin{equation}
	\begin{aligned}
		-1=&g(1) T_{2}-[r(1)+g(1)] T_{1}, \\
		\vdots & \\
		-1=&r(N) T_{N-1}+g(N) T_{N+1}-(r(N)+g(N)) T_{N}, \\
		& \vdots \\
		-1=&r\left(N_{T}\right) T_{N_{T}-1}-r\left(N_{T}\right) T_{N_{T}}.
		\label{TT}
	\end{aligned}
\end{equation}
There are two methods to obtain the mean run time, the first method is to solve the matrix directly, or we can use recursive relational equations to calculate an exact expression for the mean run time. Next we first introduce the first method.
We can rewrite  Eq.~\eqref{TT} in matrix form as
\begin{equation}
A {\bf T}=-{\bf e},
\end{equation}
where ${\bf T}=[T_{1},\ldots,T_{N},\ldots,T_{N_{T}}]^T$, ${\bf e}=[1,\ldots,1,\cdots,1]^{T}$ and 
\begin{equation}
	A=\left(\begin{array}{ccccccccc}
		-[r(1)+g(1)] & g(1) & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 \\
		\vdots & \vdots & & \vdots & \vdots & \vdots & & \vdots & \vdots \\
		0 & 0 &  \cdots& r(N) & -[r(N)+g(N)] &g(N)&\cdots&0&0 \\
		\vdots & \vdots & & \vdots & \vdots & \vdots &  & \vdots& \vdots \\
		0 & 0 & \cdots & 0 & 0 & 0 &\cdots  & r\left(N_{T}\right)&-r\left(N_{T}\right)
	\end{array}\right).
\label{matrixA}
\end{equation}
So the mean run time of myosin filament along actin is
\begin{equation}
	\left\langle T\right\rangle =\sum_{N=1}^{N_{T}} P\left( N|F, {\rm [ATP]}\right) T_{N}.
\end{equation}
In the following we present the second method of solving the run time by recursive relational equation.
From the $N$-th equation in Eq.~\eqref{TT}, we get a recursive relation
\begin{equation}
		T_{N-1}-T_{N}=\frac{g(N)}{r(N)}\left(T_{N}-T_{N+1}\right)-\frac{1}{r(N)}.
        \label{TNN}
\end{equation}
From the  last equation in Eq.~\eqref{TT}, we find
\begin{equation}
	T_{N_{T}-1}-T_{N_{T}}=\frac{-1}{r\left(N_{T}\right)}<0,
\end{equation}
Therefore, $T_{N_{T}-1}<T_{N_{T}}$.
Then according to the recursive relation Eq.~\eqref{TNN}, we can introduce
$
	T_{N-1}<T_{N},
$
Therefore the series $\left \{  T_{N}\right \} $ is monotonically increasing.
We let $A_{N}=T_{N}-T_{N-1}$, then Eq.~\eqref{TNN} above changes to
\begin{equation}
	A_{N}=\frac{g(N)}{r(N)}A_{N+1}+\frac{1}{r(N)},
\end{equation}
By further iterations, we finally get
\begin{equation}
	\begin{aligned}
		A_{N}&=\frac{g(N)}{r(N)}A_{N+1}+\frac{1}{r(N)} \\
		&=\frac{g(N)}{r(N)}\left(\frac{g(N+1)}{r(N+1)}A_{N+2}+\frac{1}{r(N+1)}\right)+\frac{1}{r(N)} \\
		&=\cdots\\
		&=\frac{g(N)g(N+1) \cdots g(N_{T}-1)}{r(N)r(N+1) \cdots r\left(N_{T}-1\right)}A_{N_{T}}+\frac{g(N)g(N+1)\cdots g\left(N_{{T}}-2\right)}{r(N)r(N+1) \cdots r\left(N_{T}-2\right)}\frac{1}{r\left(N_{T}-1\right)}+\cdots+\frac{g(N)}{r(N)}\frac{1}{r\left(N+1\right)}+\frac{1}{r(N)},
		\label{TNNN}
	\end{aligned}
\end{equation}
from the last equation of Eq.~\eqref{TT}, we know 
$A_{N_{T}}=1/r(N_{T}),$
so the above Eq.~\eqref{TNNN} can  be  simplified
\begin{equation}
		A_{N}=\sum_{j=N}^{N_{T}} \frac{1}{r(j)} \prod_{k=N}^{j-1} \frac{g(k)}{r(k)}.
\end{equation}
By summing over $A_{N}$, we can derive
\begin{equation}
	T_{N}-T_{1}=\sum_{i=2}^{N} A_{i}.
\end{equation}
In addition,
$A_{1}=T_{1}-T_{0}=T_{1},$
it is deduced that
\begin{equation}
	T_{N}=\sum_{i=1}^{N} A_{i}=\sum_{i=1}^{N}\left( \sum_{j=i}^{N_{T}} \frac{1}{r(j)} \prod_{k=i}^{j-1} \frac{g(k)}{r(k)}\right).
\end{equation}
So the mean run time of myosin filament along actin is
\begin{equation}
	\left\langle T\right\rangle =\sum_{N=1}^{N_{T}} P\left( N|F, {\rm [ATP]}\right) T_{N}=\sum_{N=1}^{N_{T}} P\left( N|F, {\rm [ATP]}\right)\sum_{i=1}^{N}\left( \sum_{j=i}^{N_{T}} \frac{1}{r(j)} \prod_{k=i}^{j-1} \frac{g(k)}{r(k)}\right).
	\label{meanT}
\end{equation}
	And the mean run length of the myosin II motors along actin is 
\begin{equation}
	\left\langle l\right\rangle =	\left\langle T\right\rangle V(F, [\rm ATP]).
	\label{meanL}
\end{equation} 
In addition, we can also calculate the probability of myosin II motors detaching  from the actin filament after moment $t$.
We can also express Eq.~\eqref{T} in matrix form. One can easily show that ${\bf F}(t)=[F_{1}(t), F_{2}(t),\ldots, F_{N}(t),\ldots, F_{N_{T}}(t)]$ satisfies 
\begin{equation}
	\frac{d {\bf F}(t)}{d t}=A {\bf F}(t)+\left(\begin{array}{c}
		r(1) \delta (t) \\
		0 \\
		\vdots\\
		0
	\end{array}\right),
\label{Ft}
\end{equation}
with $A$ given in Eq.~\eqref{matrixA}. 

Let ${\mathbf \rho}(t)=\int_{0}^{t}{\bf F}(t)dt=[\rho_{1}(t),\ldots, \rho_{N_{T}}(t)]$ represent the probability of finding the myosin filament detachment from actin
before time $t$, then from Eq.~\eqref{Ft}, we have
\begin{equation}
	{\bf F}(t)-{\bf F}(0)=A {\mathbf \rho}(t)+\left(\begin{array}{c}
		r(1)  \\
		0 \\
		\vdots\\
		0
	\end{array}\right),
	\label{rho}
\end{equation}
where ${\bf F}(0)={\bf 0}$.
Accordingly, we obtain the existence probability $\tilde{\rho}(t)={\bf e}-{\bf \rho}(t)=\int_{t}^{\infty}{\bf F}(t)dt=[\tilde{\rho}_{1}(t), \ldots, \tilde{\rho}_{N_{T}}(t)]$, which is the probability that myosin  filament detaches from actin after time $t$, with ${\bf F}(\infty)={\bf 0}$. Then we get
\begin{equation}
	\frac{d {\mathbf {\tilde \rho}}(t)}{d t}={\bf F}(\infty)-{\bf F}(t)=-{\bf F}(t).
	\label{Ftt}
\end{equation}
Thus, combining with the Eq.~\eqref{rho} and Eq.~\eqref{Ftt}, we deduce that
\begin{equation}
		-\frac{d {\mathbf {\tilde \rho}}(t)}{d t}=A\left( {\bf e}-{\bf {\tilde \rho}}(t)\right) +\left(\begin{array}{c}
		r(1)  \\
		0 \\
		\vdots\\
		0
	\end{array}\right),
\end{equation}
with $A{\bf e}^{T}=\left[ -r(1),\ldots,0\right] ^{T}$, so
\begin{equation}
	\frac{d {\mathbf {\tilde \rho}}(t)}{d t}=A{\bf {\tilde \rho}}(t),
\end{equation}
with ${\bf {\tilde \rho}}(0)={\bf e}$. 
Finally, we give the expression for existence probability 
$
	\tilde{\rho}(t)=e^{At}{\bf e}.
$
 So the mean existence probability is
\begin{equation}
	\label{meanrho_S}
	\left\langle \tilde{\rho}(t) \right\rangle =\sum_{N=1}^{N_{T}} P\left( N|F, {\rm  [ATP]}\right) \tilde{\rho}_{N}(t).
\end{equation}
Let ${\bf P}=[P\left( 1|F, {\rm  [ATP]}\right), \ldots, P\left( N|F, {\rm  [ATP]}\right),\ldots, P\left( N_{T}|F, {\rm  [ATP]}\right)]$, Eq.~\eqref{meanrho_S} can be rewritten as $	\left\langle \tilde{\rho}(t) \right\rangle={\bf P}e^{At}{\bf e}$. \\


{\centering\subsection {Sensitivity analysis and clustering results of parameters}}

\begin{figure}[htbp]
	\includegraphics[scale=0.55]{Kmeans.eps}\\
	\caption{PCA visualization of clustering results. Circular markers represent points belonging to Cluster 1, square markers represent points belonging to Cluster 2, and the two crosses indicate the Cluster Centers.}
	\label{kmeans}
\end{figure}	

The parameters we provide for theoretical predictions have a good level of robustness. When the parameters are perturbed, the overall trend of the predicted curves are not affected, and the original nature of the predicted curves remain preserved, with just slight alterations in the values.

In Eq.\eqref{omega_off_on}, the parameter $\delta$ is set to $\delta=1.3$ \cite{veigel2003load}. In the activation curve equation, Eq.\eqref{NT}, the $N_0$ is assigned a value of $N_0 = 294$ \cite{piazzesi2007skeletal}. This value represents the total number of active and inactive motors. Additionally, the $F_0$ is set to $F_0 = 480$ pN \cite{piazzesi2007skeletal}. The parameters $a_0 = 0.95$ and $b = 7.76$ in Eq.\eqref{NT} are determined by fitting the activation data from Linari $et. al.$ \cite{piazzesi2007skeletal}. The $d$ is assigned a value of $8$ nm \cite{kaya2010nonlinear}. The free energy  released by the hydrolysis of ATP  $\Delta \mu_{\rm ATP}$ is equal to $25$ $k_BT$ at [ATP] = 2000 $\mu$M \cite{wagoner2019opposing}, and is dependent on the concentration of ATP.

Currently, there is no universally accepted standard value for the worm-like chain (WLC) parameters of the neck linker in myosin II motors. The rigidity of this region can be influenced by various factors such as ion concentration, pH, and temperature. However, several studies have employed techniques like stretching experiments to measure the neck linker of myosin II motors \cite{fang2021biomechanical, li2005spectrin, fedosov2010systematic,oberhauser1998molecular}.
Combining these results, we estimate the persistence length $L_{p}$ to be approximately $0.1-1.0$ nm and the contour length $L_{c}$ to be approximately $9-40$ nm. These values serve as references, but specific values may vary depending on different experimental conditions.

In the following, we will show the fitting  procedure. 
We simultaneously fit to biophysical quantities (the velocity $V$, the average number of bound motors $\left\langle N\right\rangle $, the average force shared by each myosin $F/\left\langle N\right\rangle $, and the sliding distance $L$), which has the advantage of better constraining the fit of our model. 
Optimization of the sum of the combined errors of these four biophysical quantities to data measured in \cite{piazzesi2007skeletal}. 





\begin{table}[tbhp]
	\centering
	\caption{
		List of five sets parameter values with the smallest errors. }
	\begin{tabular}{ |p{2cm}|p{3.5cm}|p{2cm}|p{2cm}|p{2cm}|p{2cm}|p{2cm}|}
		\hline
		Parameters& Range tested &Best fit&Set 1&Set 2&Set 3&Set 4\\
		\hline
		$g_{\rm on}$ & $-10$ to 0 $k_{B}T$ \cite{wagoner2021evolution}& $-5.637$ 	&$-5.203$&$-6.938$&$-7.230$
		&$-7.151$\\
		$g_{\rm ps}$& $-20$ to $-15$ $k_{B}T$ \cite{wagoner2021evolution}& $-17.438$&$-16.240$&$-19.473$
		&$-15.721$
		&$-18.194$\\
		$g_{\rm off}$&$-10$ to 0 $k_{B}T$ \cite{wagoner2021evolution}&  $-3.499$&$-4.354$&$-2.957$
		&$-4.667$
		&$-2.947$\\
		$g_{\rm recovery}$&$-$&1.574 & 0.796&4.367&2.618
		&3.291\\
		$\omega _{\rm on}$&0.01 $-$ 10 s$^{-1}$ \cite{wagoner2021evolution}& 0.178& 0.128&0.648&0.096
		&0.189\\
		$\omega_{\rm off}$&0.1 $-$ 100000 s$^{-1}$ \cite{wagoner2021evolution}& 194.009 & 177.556&173.866&211.186
		&188.831\\
		$d^{\prime}$& 4 $-$ 8 nm \cite{kaya2017coordinated}& 4.187&4.108&4.171
		&4.355
		&4.299\\
		$L_{p}$& 0.1 $-$ 1 nm \cite{fang2021biomechanical, li2005spectrin, fedosov2010systematic,oberhauser1998molecular}&0.162&0.135&0.152&0.180
		&0.126\\
		$L_{c}$&9 $-$ 40 nm \cite{fang2021biomechanical, li2005spectrin, fedosov2010systematic,oberhauser1998molecular}& 24.561&26.628&28.010
		&21.339
		&32.611\\
		$k_{\rm pre}$& 0 $-$ 10000s$^{-1}$ \cite{kaya2017coordinated}& 5.03E$-06$ &0.000873&3.17E$-05$
		&1.50E$-05$
		&2.87E$-06$\\
		$k_{\rm mid}$& 0 $-$ 10000s$^{-1}$ \cite{kaya2017coordinated}&0.714 &1.975&0.146
		&2.582
		&0.363\\
		$\gamma$& 0.1 $-$ 1 \cite{kaya2017coordinated,kaya2010nonlinear}& 0.192&0.398&0.126
		&0.128
		&0.176\\
		\hline
	\end{tabular}\label{tab:S1}
\end{table}

We filtered a parameter dataset based on the magnitude of error, then, we selected 127 sets of parameters for further analysis.  
 Initially, we conduct a cluster analysis on the dataset, followed by visualizing the clustering results through principal component analysis. The advantage of this method is that it allows for better comprehension of the clustering results, revealing the inherent structure of the parameters.
 
We use the $k$-means clustering method to cluster parameters. This algorithm begins by randomly selecting $k$ cluster centers, each of which represents the centroid of a cluster. Then, for each data point in the dataset, the algorithm assigns it to the cluster with the nearest centroid based on a distance metric such as Euclidean distance. Once all the data points have been assigned to a cluster, the algorithm recalculates the centroids of each cluster based on the mean of all the points assigned to that cluster. This process is repeated iteratively until the centroids no longer change. After clustering the dataset into k clusters using $k$-means clustering, we obtain the index of each data point's corresponding cluster and the cluster centers.

Selecting different initial cluster centers leads to divergent clustering results. In order to select a suitable initial center, we first use hierarchical clustering to cluster the parameter dataset, terminating at $k$ clusters. Then, we select the point closest to the center from each cluster as the initial cluster center for $k$-means clustering.

The number of clusters $k$ in $k$-means clustering needs to be specified beforehand, but in practical applications, the optimal value of $k$ is often unknown. One solution to this problem is to try clustering with different values of $k$, evaluate the quality of the resulting clusters, and infer the optimal $k$ value. The quality of the clustering results can be measured by the average diameter of the clusters. 
For each cluster, we find all the data points within it, calculate the distance matrix between them, and find the maximum value in the distance matrix, which represents the diameter of the cluster. The average of all cluster diameters is the average diameter.

Generally, as the number of clusters decreases, the average diameter of the clusters increases, and when the number of clusters exceeds a certain value, the average diameter of the clusters does not change. This value of $k$ is precisely the optimal value. Thus, our optimal value of  $k$ is $k=2$, obtained using this method.

After completing the clustering process, we perform principal component analysis (PCA) on the dataset and visualize the results. PCA visualization helps us better understand the clustering outcomes. By using the PCA algorithm, we transform the clustering results into principal components and project the original parameter onto these components to obtain a new set of low-dimensional features, which contain the most informative parts of the data.

We create a scatter plot in which each point represents a sample, with its position determined by its score on the first two principal components. We utilize different colors and shapes to represent different clusters, thereby providing a clear and concise representation of the clustering outcomes, as shown in Fig.~\ref{kmeans}.


 From all the parameter sets, we select the five parameter sets, the values of these parameter sets are given in Tab.~\ref{tab:S1}. 
The parameter sets 1 and 2 belong to the Cluster 1, while the parameter sets 3 and 4, as well as the best-fit set, belong to the Cluster 2. We plotted all the theoretical prediction curves for these five parameter sets. 
By comparing the prediction curves of these five parameter sets, we find that
the perturbation of parameters has a limited influence on the trend of predicted values.
The curves predicted with these five parameter sets are similar, and each set of parameter sets corresponds to the same curve properties, with only differences in values. The tendency of curves may appear sooner or later, but they will always appear. As a result, our theoretical predictions are very reliable.
Finally, we selected the best-fit set, which resulted in the best fitting effect. The final parameters fit curve resembles the data quite closely, see Fig.~3.


Due to the discontinuity of $N_{T}$ caused by our use of integer values, the mean run time $\left\langle t \right\rangle $ values predicted with the parameters provided in the main text may be slightly inaccurate compared to the exact values; nevertheless, the trends of increasing run time will definitely occur. We have acknowledged this issue and conducted multiple extensive searches to validate our results. The figures presented herein are intended to serve as a reference for readers.

%\bibliographystyle{unsrt}
%\bibliography{mybibtex2023}

	\begin{thebibliography}{48}
	\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
	\expandafter\ifx\csname bibnamefont\endcsname\relax
	\def\bibnamefont#1{#1}\fi
	\expandafter\ifx\csname bibfnamefont\endcsname\relax
	\def\bibfnamefont#1{#1}\fi
	\expandafter\ifx\csname citenamefont\endcsname\relax
	\def\citenamefont#1{#1}\fi
	\expandafter\ifx\csname url\endcsname\relax
	\def\url#1{\texttt{#1}}\fi
	\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
	\providecommand{\bibinfo}[2]{#2}
	\providecommand{\eprint}[2][]{\url{#2}}

	\bibitem[{\citenamefont{kaya}(2010)}]{kaya2010nonlinear}
\bibinfo{author}{
	\bibfnamefont{M.} \bibnamefont{Kaya} \bibnamefont{and}
	\bibfnamefont{H.} \bibnamefont{Higuchi}},
\bibinfo{journal}{Science} \textbf{\bibinfo{volume}{329}}, \bibinfo{pages}{686}
(\bibinfo{year}{2010}).

		\bibitem[{\citenamefont{fang}(2021}]{fang2021biomechanical}
\bibinfo{author}{
	\bibfnamefont{Y.} \bibnamefont{Fang},
	\bibfnamefont{Y.} \bibnamefont{Hu},
	\bibfnamefont{F.} \bibnamefont{Cheng},
	\bibnamefont{and}
	\bibfnamefont{Y.} \bibnamefont{Xin}},
\bibinfo{journal}{Physical Review E} \textbf{\bibinfo{volume}{104}}, \bibinfo{pages}{064403}
(\bibinfo{year}{2021}).

		\bibitem[{\citenamefont{fedosov}(2010}]{fedosov2010systematic}
\bibinfo{author}{
	\bibfnamefont{D.~A.} \bibnamefont{Fedosov},
	\bibfnamefont{B.} \bibnamefont{Caswell},
	\bibnamefont{and}
	\bibfnamefont{G.~E.} \bibnamefont{Karniadakis}},
\bibinfo{journal}{Computer Methods in Applied Mechanics and Engineering} \textbf{\bibinfo{volume}{199}}, \bibinfo{pages}{1937}
(\bibinfo{year}{2010}).

		\bibitem[{\citenamefont{li}(2005}]{li2005spectrin}
\bibinfo{author}{
	\bibfnamefont{J.} \bibnamefont{Li},
	\bibfnamefont{M.} \bibnamefont{Dao},
	\bibfnamefont{C.~T.} \bibnamefont{Lim},
	\bibnamefont{and}
	\bibfnamefont{S.} \bibnamefont{Sureshs}},
\bibinfo{journal}{Biophysical Journal} \textbf{\bibinfo{volume}{88}}, \bibinfo{pages}{3707}
(\bibinfo{year}{2005}).

		\bibitem[{\citenamefont{piazzesi}(2007)}]{piazzesi2007skeletal}
\bibinfo{author}{
	\bibfnamefont{G.} \bibnamefont{Piazzesi},
	\bibfnamefont{M.} \bibnamefont{Reconditi},
	\bibfnamefont{M.} \bibnamefont{Linari},
	\bibfnamefont{L.} \bibnamefont{Lucii},
	\bibfnamefont{P.} \bibnamefont{Bianco},
	\bibfnamefont{E.} \bibnamefont{Brunello},
	\bibfnamefont{V.} \bibnamefont{Decostre},
	\bibfnamefont{A.} \bibnamefont{Stewart},
	\bibfnamefont{D.~B.} \bibnamefont{Gore},
	\bibfnamefont{T.~C.} \bibnamefont{Irving},
	\bibfnamefont{M.} \bibnamefont{Irving},
	\bibnamefont{and}
	\bibfnamefont{V.} \bibnamefont{Lombardi}},
\bibinfo{journal}{Cell} \textbf{\bibinfo{volume}{131}}, \bibinfo{pages}{784}
(\bibinfo{year}{2007}).

	\bibitem[{\citenamefont{kaya}(2017)}]{kaya2017coordinated}
\bibinfo{author}{
	\bibfnamefont{M.} \bibnamefont{Kaya},
	\bibfnamefont{Y.} \bibnamefont{Tani},
	\bibfnamefont{T.} \bibnamefont{Washio}, 
	\bibfnamefont{T.} \bibnamefont{Hisada}, \bibnamefont{and}
	\bibfnamefont{H.} \bibnamefont{Higuchi}},
\bibinfo{journal}{Nature Communications} \textbf{\bibinfo{volume}{8}}, \bibinfo{pages}{1}
(\bibinfo{year}{2017}).

		\bibitem[{\citenamefont{erdmann}(2012}]{erdmann2012stochastic}
\bibinfo{author}{
	\bibfnamefont{T.} \bibnamefont{Erdmann}
	\bibnamefont{and}
	\bibfnamefont{U.~S.} \bibnamefont{Schwarz}},
\bibinfo{journal}{Physical review letters} \textbf{\bibinfo{volume}{108}}, \bibinfo{pages}{188101}
(\bibinfo{year}{2012}).

		\bibitem[{\citenamefont{wagoner}(2021)}]{wagoner2021evolution}
\bibinfo{author}{
	\bibfnamefont{J.~A.} \bibnamefont{Wagoner} \bibnamefont{and}
	\bibfnamefont{K.~A.} \bibnamefont{Dill}},
\bibinfo{journal}{Proceedings of the National Academy of Sciences} \textbf{\bibinfo{volume}{118}}, \bibinfo{pages}{e2101871118}
(\bibinfo{year}{2021}).


	\bibitem[{\citenamefont{veigel}(2003}]{veigel2003load}
\bibinfo{author}{
	\bibfnamefont{C.} \bibnamefont{Veigel},
	\bibfnamefont{J.~E.} \bibnamefont{Molloy},
	\bibfnamefont{S.} \bibnamefont{Schmitz},
	\bibnamefont{and}
	\bibfnamefont{J.} \bibnamefont{Kendrick-Jones}},
\bibinfo{journal}{Nature cell biology} \textbf{\bibinfo{volume}{5}}, \bibinfo{pages}{980}
(\bibinfo{year}{2003}).

	\bibitem[{\citenamefont{wagoner}(2019}]{wagoner2019opposing}
\bibinfo{author}{
	\bibfnamefont{J.~A.} \bibnamefont{Wagoner}
	\bibnamefont{and}
	\bibfnamefont{K.~A.} \bibnamefont{Dill}},
	\bibinfo{journal}{Molecular biology and evolution} \textbf{\bibinfo{volume}{36}}, \bibinfo{pages}{2813}
	(\bibinfo{year}{2019}).

	\bibitem[{\citenamefont{oberhauser}(1998}]{oberhauser1998molecular}
\bibinfo{author}{
	\bibfnamefont{A.~F.} \bibnamefont{Oberhauser},
	\bibfnamefont{P.~E.} \bibnamefont{Marszalek},
	\bibfnamefont{H.~P.} \bibnamefont{Erickson},
	\bibnamefont{and}
	\bibfnamefont{J.~M.} \bibnamefont{Fernandez}},
\bibinfo{journal}{Nature} \textbf{\bibinfo{volume}{393}}, \bibinfo{pages}{181}
(\bibinfo{year}{1998}).


	\end{thebibliography}





	















\end{widetext}
\end{document}
