Distribution of Myosin Attachment Times Predicted from Viscoelastic Mechanics of Striated Muscle

We demonstrate that viscoelastic mechanics of striated muscle, measured as elastic and viscous moduli, emerge directly from the myosin crossbridge attachment time, tatt, also called time-on. The distribution of tatt was modeled using a gamma distribution with shape parameter, p, and scale parameter, β. At 5 mM MgATP, β was similar between mouse α-MyHC (16.0 ± 3.7 ms) and β-MyHC (17.9 ± 2.0 ms), and p was higher (P < 0.05) for β-MyHC (5.6 ± 0.4 no units) compared to α-MyHC (3.2 ± 0.9). At 1 mM MgATP, p approached a value of 10 in both isoforms, but β rose only in the β-MyHC (34.8 ± 5.8 ms). The estimated mean tatt (i.e., pβ product) was longer in the β-MyHC compared to α-MyHC, and became prolonged in both isoforms as MgATP was reduced as expected. The application of our viscoelastic model to these isoforms and varying MgATP conditions suggest that tatt is better modeled as a gamma distribution due to its representing multiple temporal events occurring within tatt compared to a single exponential distribution which assumes only one temporal event within tatt.


Introduction
Advances in optical techniques have allowed detailed analysis of the isolated single myosin crossbridge. Upon formation of a myosin crossbridge, isomerization of the myosin molecule provides a unitary force proportional to the length of the lever arm and the extent of the lever arm swing, also called the power stroke [1][2][3][4]. The time duration of crossbridge attachment, sometimes called time-on, is also measured by optical trapping. The distribution of attachment times has been found to depend upon the myosin isoform, nucleotide availability, point mutations, and several other factors [5][6][7][8][9][10][11][12]. This time duration of crossbridge attachment plays a significant role in determining muscle performance observed as force and velocity of contraction at the single fiber level.
There remains, however, a considerable challenge to understand and describe how the unitary force and the temporary attachment times of myosin molecules are manifested at the level of muscle tissue, which possesses a threedimensional lattice structure and can be studied intact or submerged in physiological solutions after removal of the plasma membrane. Viscoelastic mechanics of muscle tissue represents one such macroscopic consequence of these molecular phenomena. We describe in this paper a quantitative justification and methodology for estimating the distribution of myosin crossbridge attachment times based on viscoelastic mechanics measured in striated muscle fibers.
Length perturbation analysis of muscle provides a means to quantify viscoelastic mechanics and entails applying a small length change at one end of a muscle and recording the force response on the other end. The length perturbation and force response are then subjected to Fourier transformation, and the dynamic mechanical properties of the muscle are characterized by the complex ratio of the transformed force response normalized to cross-sectional area (or σ(ω) = Fourier transformed tensile stress) divided by the transformed length perturbation normalized to initial length (or ε(ω) = Fourier transformed strain). This complex ratio is also called the mechanical transfer function or complex modulus of the muscle, Y(ω) = σ(ω)/ε(ω). The real and imaginary parts are, respectively, termed the elastic modulus and viscous modulus of the muscle [13,14]. The frequency characteristic dips (low values) and shoulders (high values) observed in the elastic and viscous moduli of activated muscle are sensitive to species of origin, myosin isoform, and varying concentrations of MgATP and therefore reflect the macroscopic consequence of myosin crossbridge kinetics [4,11,[15][16][17][18]. The dips in the moduli, most notably the negatively valued viscous modulus, occur near the frequencies at which the muscle operates in vivo [14]. The shoulders of the moduli, however, are more prominent in magnitude. According to our previous modeling endeavor [19], the shoulders appearing at these higher frequencies of the moduli reflect the mechanical consequences of intermittent myosin crossbridge formation. In that work, we proposed that a two-state model of the acto-myosin crossbridge governed by first-order kinetics gives rise to a viscoelastic work-absorbing property, termed the C-process by Kawai and colleagues [13,16,20], which is characterized by an exponential rate constant equivalent to the myosin crossbridge off rate termed g by Huxley [21]. The mean myosin attachment time, based on a single exponential distribution of attachment times, could be estimated as the reciprocal of this exponential rate constant, 2πc, after fitting (1) to a measured complex modulus [14,19]: While the single exponential representation of C-process and its interpretation have been valuable for examining a variety of muscle types under a number of conditions [7,11,13,14,19,20,22], the assumptions of first-order kinetics and a single exponential distribution of myosin attachment times are limiting. It is known, for example, that multiple time periods associated with multiple biochemical states make up the myosin crossbridge cycle. These multiple states constitute the entire myosin crossbridge attachment time. Such an addition of multiple smaller time periods would result in a multiple exponential distribution or a gamma distribution of the crossbridge attachment times rather than a single exponential distribution assumed by first-order kinetics.
In the present study, we examine the mechanical consequences of the distribution of myosin attachment times represented by a gamma distribution, which allows the description of more discrete and longer lived time-on and yet also allows for the possibility of a single exponential function. The consideration of the gamma distribution effectively poses the hypothesis that multiple temporal events which occur within the time of attachment can be discerned in the viscoelastic response of striated muscle at the macroscopic level. We provide an analytical solution to the mechanical consequences that emerge as the shoulders of the elastic and viscous moduli. We also demonstrate the validity of the analytical solution with computer simulations. Finally, we demonstrate the application of the gamma distribution representation of the C-process in a comparison of mouse cardiac myosin heavy chain isoforms, α-MyHC and β-MyHC, subjected to varying concentrations of MgATP.

Mathematical Modeling
Our goal in this section is to derive a mathematical expression that describes the viscoelastic component of a force response to a length perturbation that has been applied to an ensemble of myosin molecules intermittently attaching and detaching from actin within a half sarcomere of striated muscle. The force response and length perturbation will form the bases for a mechanical transfer function, which is equivalent to a viscoelastic complex modulus measured in frequency space. We use a two-state model of myosin attachment and detachment as our starting point in modeling the mechanical consequences of intermittently attached actomyosin crossbridges. It is important to point out that in the following development we do not assume first-order kinetics as governing myosin attachment and detachment; in doing so we will be able to provide a more generalized model of the viscoelastic mechanics due to the myosin crossbridges.
The time periods of attachment, t att , and detachment, t det , are considered here to be random variables whose values are governed by stochastic processes independent of force, stress, strain, length, velocity, and each other. The total time period for a myosin crossbridge cycle is also a random variable, t cycle = t att + t det . All myosin heads are assumed independent of each other. During detachment, force due to an individual myosin head is zero. During attachment, force is given as the product of the unitary force due to the power stroke, the length displacement of an elastic element in series with the crossbridge, and the stiffness of that elastic element ( Figure 1). We assume that length perturbations are sufficiently small to elicit a linear force response of the elastic element in series with the crossbridge.
The elastic element represents the most compliant structures between the M-and Z-lines of the sarcomere in series with the attached crossbridge. Previous work by others has identified the most compliant structures as the head and neck regions of the myosin S1 segment including the lever arm bound with the essential and regulatory light chains and that portion of myosin S2 segment not incorporated into the thick filament backbone [4,23,24]. We also assume that when the myosin is in a postpower stroke state the most compliant elements are taut, that is, not slack as might occur if the myosin bound to actin at a point too close to the M-line to stretch the S2 segment [4,23,24].

Viscoelastic Mechanical Response of Half Sarcomere.
Here we provide a theoretical description of the force response of a half sarcomere to an externally applied length perturbation. We assume here that the total force that might be recorded from a half sarcomere, F total (t), is given as the sum of the forces produced by each attached myosin crossbridge due to its unitary force, F uni , and its dynamic response due to the length perturbation, f (t) (Figure 1). The total force produced by a summation of attached crossbridges will be denoted as Journal of Biomedicine and Biotechnology where M = the number of attached crossbridges in the half sarcomere at time t. When an external length perturbation is applied, the dynamic response of the ith attached crossbridge is given as the stiffness-displacement product, which is a linear approximation suitable for very small length perturbations: where k stiff = the stiffness coefficient of the crossbridge elastic element, L hs (t) = displacement imposed upon of the half sarcomere by an externally driven perturbation, and t ini = the instant of initial crossbridge attachment. Under isometric conditions, the total force is simply equal to the total number of attached crossbridges, M, multiplied by the mean unitary force, F uni . The dynamic force recorded over the half sarcomere, F hs (t), can then be defined as the total force minus the isometric force and written as follows: We now use the Inverse Fourier Transform definitions for F hs (t) and L hs (t), namely, where F hs (ω) and L hs (ω) are the Fourier Transform representations of F hs (t) and L hs (t), respectively; then (4) can be written as follows: We now define the variable, τ = t -t ini , as the time period between time, t, and the instant of the most recent attachment of a crossbridge, t ini . For the ith crossbridge, t ini(i) = t -τ (i) . The term (e iωt − e iωtlni(i) ) in (6) then becomes (1 − e −iωτ(i) )e iωt and we can remove the Inverse Fourier Integrals of (6): The summation over a large number of attached crossbridges can be replaced as the product of M and the expected value of the summed terms. To do so would require our providing the probability density functions (PDFs) for those random variables represented as bearing the i-subscript in (7). We assume here that the random variables k stiff and τ are independent of each other. Thus, the result for the expected value of k stiff would be independent of its probability density function and equal to the mean value for k stiff . Equation (7) can be written as where k stiff = mean stiffness of the elastic element and PDF τ (t) = the probability density function for the random variable τ as a function of time, t. Equation (8) represents the force response in frequency space due to an externally applied length perturbation on striated muscle with no strain dependence on myosin crossbridge kinetics. The number of crossbridges attached at any time, M, can be replaced by the total number of myosin heads available to form crossbridges, N, multiplied by the probability that any one myosin head has formed a crossbridge, that is, the duty ratio t att /t cycle . The mechanical transfer function, which represents a measured viscoelastic complex modulus, can then be written as follows: where t att = mean time duration of crossbridge attachment and t cycle = mean time of a complete crossbridge cycle. Equation (9) represents a generalized mathematical description of the visco-elastic response of intermittently attaching myosin crossbridges. The form of (9) suggests that the force response is dictated by PDF τ (t), so long as the other random variables are sufficiently represented by their respective means. A physical interpretation of PDF τ (t) and its relation to the probability density function for t att , PDF tatt (t), are presented below.

Survival Function.
As defined above, the random variable τ represents the time period between any given time t and the instant any bound crossbridge was formed (Figure 1(b)). In other words, τ represents the time period any bound crossbridge has survived up to time t. The probability (Pr) that a crossbridge has survived for a time period τ, we will call it S τ (t), is defined as the probability that τ is less than the time of crossbridge attachment, t att [25]. The probability density function of crossbridge attachment time PDF tatt (t) and the survival function S τ (t) are related as follows: The PDF τ (t) needed in (8) is simply a normalization of S τ (t) and is defined as the survival function of (10a) divided by its integral: With (9) we have provided a mathematical representation of the energy-absorbing viscoelastic complex modulus that arises from intermittently attached myosin crossbridges, and we have not assumed a specific scheme for the biochemical steps that govern the distributions of the random variable t att . Thus, any measure of the viscoelastic complex modulus, given by the left hand side of (9), can be used to calculate an estimate of PDF τ (t) and by extension through (10b) also be used to calculate an estimate of PDF tatt (t).

Gamma Distribution Representation of PDF tatt (t).
In our previous work we chose PDF tatt (t) to be represented by a single exponential distribution [19], which by extension of (10a), (10b), and (11) also defined PDF τ (t) as the same single exponential distribution. That choice reflected an assumption of first-order kinetics governing the lifetime of the myosin crossbridge in the attached state. The result was a mathematical representation of the mechanical consequences of length perturbation that was equivalent to the historically represented C-process shown in (1).
Here we do not assume first-order kinetics. Instead, we consider PDF tatt (t) to be represented by a gamma distribution, which is used to describe random variables that represent a time period made up of the sum of several smaller time periods [26,27] where p = a shape parameter reflecting the number of smaller time periods summed together to produce total time attached, t att , β = a scale parameter reflecting the average duration of the smaller time periods. The mean t att is calculated as pβ, and Γ(p) = the gamma function evaluated at p, a normalization factor. The two parameters p and β are sufficient to provide a gamma distribution description of PDF tatt (t). The case of p = 1 is a special case that is equivalent to a single exponential function with mean β. Several examples of gamma distributions, including a case of p = 1, are illustrated in Figures 2(a) and 2(b), which demonstrate the effectiveness of the gamma distribution representation in accommodating a wide variety of probability density functions for t att . Upon applying (10a), we attain the associated survival function S τ (t): where Γ(p, t/β) = the upper incomplete gamma function [26].
Upon applying (13a) and (13b) to (11) we have We can now state that the complex modulus that would arise from an ensemble of myosin crossbridges, whose distribution of t att can be represented by a gamma distribution presented in (12), is given by the insertion of PDF τ (t) of (14) into the expression for the mechanical transfer function of Time (ms) The two parameters of the gamma distribution (p and β) permit description of a wide variety of distributions useful for representing the probability density function for myosin crossbridge attachment time, PDF tatt (t). Notably, a single exponential results when p = 1 (a), and a Gaussian distribution is approximated as p increases.
(c, d) The survival function, S τ (t), refers to the probability that an attached crossbridge will survive to time t. The analytical relationship between the PDF tatt (t) and S τ (t) is provided by (10a) and (10b). (9). Details of the evaluation of the integral in (9) after such a substitution are provided in the Appendix. The result is where ξ = N(t att /t cycle )k stiff . Equation (15) describes the predicted complex modulus that would arise from length perturbation analysis of an ensemble of myosin crossbridges whose distribution of attachment times is represented by a gamma distribution. From a measured viscoelastic complex modulus represented on the left-hand side of (15) we could use nonlinear leastsquares methods to estimate the values of the three independent parameters ξ, p, and β, on the right-hand side of (15).

Computer Simulations.
The force response of a virtual half sarcomere was simulated for a time period of two seconds for many (20,000) independent myosin heads alternately attaching and detaching to actin according to independent stochastic processes governing the random variables t att and t det . For any one computer simulation, random numbers representing t att were generated to conform 6 Journal of Biomedicine and Biotechnology with a gamma distribution with values for p ranging from 1 to 48 and for β ranging from 2 to 24 such that the mean time attached was either 24 ms or 96 ms. Random numbers representing t det were generated with values for p = 2 and β = 48 ms such that the mean time detached was always 96 ms. Figure 3(a) illustrates the PDF tatt (t) when p = 3 and β = 8. Each random number generated from this PDF tatt (t) was used to dictate the time period over which a virtual crossbridge was attached for each occasion, as indicated in Figure 3(b). Figure 3(c) illustrates the PDF for t det , PDF tdet (t), with p = 2 and β = 48. Each random number generated from this PDF tdet (t) was used to dictate the time over which a virtual crossbridge was detached for each occasion. Figure 3(b) illustrates an example of successive attachment and detachment time periods.
Sinusoidal length perturbations of the virtual half sarcomere were simulated as having amplitude 1 nm and frequencies over the range 1-250 Hz. Figure 3(d) illustrates an example perturbation at 1 Hz over a one-second time interval. Each crossbridge was assigned a stiffness constant of 1 pN/nm [28]. The change in length of the crossbridge relative to the time of initial attachment was multiplied by the stiffness constant to simulate the force resulting from the strain on the elastic element of the crossbridge. Figure 3(e) illustrates an example series of successive attachments and detachments over one second. Figure 3(f) illustrates the force deflection that would occur for a crossbridge during each time of attachment. The resulting force deflections of 20,000, independent crossbridges were summed to provide an equivalent two-second period for the ensemble in the virtual half sarcomere. An example of the resultant force from 20,000 crossbridges is given in Figure 3(g). This simulated force of the half sarcomere, F hs (t), was fit using a simplex method to a sine function, whose amplitude (Amp) and phase (φ) permitted the calculation of two components that were in phase and out of phase with respect to the length perturbation. The value of each moduli was then calculated from the amplitude and phase as elastic modulus = Amp cos(φ) and viscous modulus = Amp sin(φ). The smooth line in Figure 3(g) presents an example fit of a sine function to the resultant force, and Figure 3(h) presents the in-phase and out-of-phase components provided by the Simplex method. The resulting elastic and viscous moduli, given in units of pN/nm due to the lack of normalization factors used in the simulation, were fit simultaneously to (15) using a nonweighted Levenberg-Marquardt nonlinear leastsquares routine. The diagonal of the resulting covariance matrix indicated the variances of parameter estimates. The square root of the variance thus provided the standard error for each parameter estimate, which was then multiplied 1.96 to provide the 95% confidence range. All computer simulations, random number generation, curve-fittings, and parameter estimation routines were performed using IDL Version 7.0 (ITT, Boulder, CO).

Solutions.
All reagents were purchased from Sigma (St. Louis, Mo). Solutions were formulated by solving equations describing ionic equilibria [29]. Concentrations are expressed in mmol/L unless otherwise noted. Relaxing solu- Mouse left ventricular skinned myocardial strips were prepared using methods similar to those described previously to yield thin strips (∼140 μm diameter, ∼800 μm length) with longitudinally oriented parallel fibers [7]. These strips were chemically skinned for 2 hr at 22 • C and stored at −20 • C for no more than 5 days. At the time of study, aluminum T-clips were attached to the ends of a strip ∼150 μm apart. The strip was mounted between a piezoelectric motor (Physik Instrumente, Auburn, MA) and a strain gauge (SensoNor, Horten, Norway), lowered into a 30 μL droplet of relaxing solution maintained at 37 • C and incrementally stretched to and maintained at 2.2 μm sarcomere length detected by videography and digital Fourier transform techniques (IonOptix, Milton, MA).
Strips were calcium activated at pCa 4.5 and subjected to decreasing concentrations of 5, 2, and 1 mM MgATP by exchanging equal volumes of rigor solution. Sinusoidal perturbations of amplitude 0.125% strip length were applied over the frequency range 0.125-250 Hz. The elastic and viscous moduli were calculated from the recorded tension transient as the relative magnitudes of the in-phase and out-of-phase components with respect to the imposed sinusoidal length perturbations [13,20,30]. The measured complex modulus was fit to (1), representing the single exponential distribution model of t att , and to (16) below, representing the gamma distribution model of t att , using a non weighted Levenburg-Marquardt non linear least-squares routine running within IDL Version 7.0 (ITT, Boulder, CO) Mean myosin time attached was calculated as (2πc) −1 when the fit was performed with (1) and as pβ for fits using (16). The correlation coefficients between the recorded data and the fitted models were calculated as a Pearson correlation coefficient. The two models used to calculate the mean myosin time attached were compared using linear correlation. Data points are presented as mean ± sem.

Computer Simulations.
The computer simulations presented here served two purposes: first, to demonstrate the accuracy of the analytical derivation resulting in (15) and second, to demonstrate that estimating the parameters of (15) by the nonlinear least-squares methods provided reasonable estimates of the parameters known a priori to underlie the computer generated elastic and viscous moduli.  (15) to the simulated data provided a comparison between the PDF tatt (t) used to generate the data and the gamma distribution predicted from the fits. In this case, the predicted gamma distributions are nearly indistinguishable from the original PDF tatt (t).
(d) The survival function, S τ (t), was also predicted from the parameter estimates and again very closely resembles the original S τ (t). and β = 3, and p = 24 and β = 4. The thin solid line represents the best fit of (15) to the simulated data. It is interesting to note that the first example, p = 1 and β = 24, represents the results from a single exponential model of PDF tatt (t). The elastic modulus monotonically rises to a maximum, and the viscous modulus monotonically rises to a max at 6.6 Hz and then monotonically falls. It is important to note that the shape of the viscous modulus when p = 1 is symmetrical about the peak when shown on the log-scaled axis. compared to the single exponential distribution. Also, the viscous modulus is no longer symmetrical about its peak. Table 1 provides a comparison of parameter values used to generate the simulations and the parameter estimates calculated by the nonlinear least-squares routine fitting (15). For the three fitted parameters p, β, and ξ, the 95% confidence intervals for the parameter estimates were narrow compared to the magnitude of the values. These confidence intervals demonstrate the high precision and uniqueness of the parameter estimates. In addition, the correlation coefficients between values used and values estimated were consistently greater than 0.997. Based on the narrow confidence intervals for parameter estimates and the high correlation Table 1: Parameter values used to generate random numbers in computer simulations of a force response due to myosin crossbridges subjected to sinusoidal length perturbations, and estimated parameter values from that force response fit to (15). For all simulations t det = 96 ms, N = 20, 000 and k stiff = 1 pN/nm. Est[] = estimated value of bracketed parameter, () indicates the ±95% confidence intervals for estimated value. The correlation between actual and estimated values resulted in a correlation coefficient greater than 0.997 for each parameter. between parameter estimates and actual parameter values, we are confident that the analytical expression provided by (15) accurately represents the mechanical consequences of temporarily attached myosin crossbridges, whose times of attachment are gamma distributed.
For the three examples given in Figures 4(a) and 4(b), the estimated parameter values were used to calculate the PDF tatt (t) and S τ (t) using (12) and (13a), respectively. Figure 4(c) illustrates a comparison between the PDF tatt (t) used to generate the simulation and that calculated from the parameter estimates. Figure 4(d) provides a similar comparison for S(τ). For both PDF tatt (t) and S τ (t), the estimated parameters reasonably reproduced the original functions used to generate the simulations. These comparisons illustrate the robustness with which parameters of the gamma distribution can be estimated using the nonlinear least-squares methods, which fit (15) to the elastic and viscous moduli.

Viscoelastic Mechanics.
Muscle strips from both groups activated to a similar maximum developed tension of 20.7 ± 1.8 mN·mm −2 for WT (n = 6) and 22.1 ± 1.7 mN·mm −2 for PTU (n = 6). Figure 5 illustrates the elastic and viscous moduli recorded for maximally activated muscle strips isolated from WT and PTU mice, representing, respectively, the α-MyHC (Figures 5(a) and 5(b)) and β-MyHC (Figures 5(c) and 5(d)) isoforms. The myosin crossbridge kinetics in the α-MyHC are faster than those in β-MyHC and are reflected in the higher frequencies ranges for the dips and shoulders in the α-MyHC. For example, the most prominent dip and shoulder of the viscous modulus at 5 mM MgATP occur at ∼9 and 90 Hz, respectively, in the α-MyHC and occur at lower frequencies 2 and 30 Hz in the β-MyHC. These observations are consistent with a longer-lived t att in the β-MyHC compared to α-MyHC [5,9].
As MgATP is reduced the dips and shoulders of the moduli shift to lower frequencies, as best seen in the viscous modulus. The muscle also becomes stiffer, which is reflected in higher values of the elastic modulus, due to a greater fraction of crossbridges formed at any one time. These observations with lowering MgATP concentrations are consistent with a prolonged t att due to a prolonged rigor state prior to MgATP binding to myosin and subsequent myosin detachment from actin. (16) to the elastic and viscous moduli provided estimates of the gamma distribution parameters, p and β, and a calculation of mean t att as the product pβ. At 5 mM MgATP the shape parameter p was lower in the α-MyHC compared to β-MyHC (Figure 6(a)). At lower MgATP (1 and 2 mM), the shape parameter rises in both myosin isoforms and is not different between them. At all MgATP conditions and in both isoforms, the estimate for p was consistently near or greater than 2. This finding suggests that a single exponential distribution, which would have been represented by the gamma distribution with p = 1, did not fit the data as well as a gamma distribution with p > 2. Our finding p > 2 suggests two or more intermediate time periods within t att . These results imply that the PDF tatt (t) reflected in the recorded elastic and viscous moduli is better represented by a distribution that describes multiple time intervals within t att . In the strictest interpretation of the gamma distribution shape parameter, the rise in p with decreasing MgATP implies additional elemental time periods summed together to produce the total t att . We would, however, caution that p approaching 10 as MgATP is lowered may be more reflective of an increasing t att and not suggestive of additional biochemical states, as the lower MgATP would prolong only one biochemical state, the rigor state [4].

Parameters of Gamma Distribution Model of t att . Fits of
The parameter β effectively represents the duration of the intermediate time period that is summed p times to make up the total t att . As illustrated in Figure 6(b), parameter β was found to be on the order of 1.5-2 ms for both α-MyHC and β-MyHC, except at 5 mM MgATP where β was 3.5 ms in the β-MyHC. We expected this parameter to have increased with decreasing MgATP, as lowering MgATP prolongs the rigor state. This parameter, however, was not as sensitive as p in describing the prolongation of t att with decreasing MgATP. Fluctuations and asymmetries in the elastic and viscous moduli, which are characteristic of the gamma distribution for PDF tatt (t), may not be easily seen in the averaged values for the moduli presented in Figure 5. To better demonstrate these fluctuations, an example pair of elastic and viscous moduli recorded from one β-MyHC muscle preparation at 1 mM MgATP is presented in Figures 6(c) and 6(d). The A-and B-processes have also been subtracted from the measured elastic and viscous moduli resulting in the Cprocess, which reflects the viscoelastic mechanics modeled here. The dotted lines shown in Figures 6(c) and 6(d) represent the expected C-process when p = 1, the single exponential distribution. The recorded C-process, however, demonstrates fluctuations in the elastic modulus and an asymmetry in the viscous modulus different from the p = 1 case and consistent with the fluctuations and asymmetries for p > 1 illustrated in Figure 4.
The mean t att was prolonged in the β-MyHC compared to α-MyHC at each MgATP examined and was also prolonged as MgATP was reduced (Figure 7(a)). These findings are consistent with known consequences of myosin isoform and MgATP availability [5,9,16]. The calculation of the mean t att based on the gamma distributed t att was generally higher than that calculated based on the exponential distributed t att , showing that the shape of the PDF tatt (t) affects the calculated value of mean t att . Notably, the correlation between t att based on the gamma distribution and that based on the exponential * (c) (d) Figure 6: Results of fitting elastic and viscous modulus to gamma distributed and single exponential distributed models of PDF tatt (t). (a, b) The shape parameter, p, and the scale parameter, β, were estimated from fitting the gamma distribution model of (16) to the measured elastic and viscous moduli. * P < 0.05 by t-test between MyHC isoforms at same MgATP. (c) The elastic modulus recorded from one PTU β-MyHC muscle preparation at 1 mM MgATP is shown with its fit to (16), that is, the gamma distribution (GD) model. The goodness of fit is representative of the fits for all samples in this study. The A-and B-processes of (16) were then subtracted from the recorded elastic modulus to demonstrate fluctuations in the remaining C-process, which are more apparent when compared against the single exponential model (p = 1) of the t att distribution (dotted line). These fluctuations are mechanical consequences of the gamma distributed t att . (d) The corresponding viscous modulus was also subjected to fit and subtraction of A-and B-processes. The resulting C-process is compared against the single exponential model of the t att distribution (dotted line). The asymmetry in the peak of the viscous modulus is subtle, but apparent at frequencies higher than 11 Hz. distribution model was very strong, r 2 = 0.932 (Figure 7(b)), which indicates that both methods are capable of making comparisons of t att among multiple groups or conditions. The Pearson correlation coefficient between the recorded data and the fitted models was very high, but was higher in the gamma distribution model (r 2 = 0.986 ± 0.002) compared to single exponential model (r 2 = 0.980 ± 0.002). We would expect the estimate of mean t att based on the gamma distribution to be somewhat more accurate than those based on the single exponential distribution. Unfortunately, at this time and without some other independent method to verify the actual mean t att in the muscle strip, we cannot say with certainty what bias is introduced into our estimate of t att using either model. The parameter estimates resulting from the fits to (16) were used to predict the distributions of t att , that is, PDF tatt (t), for each MyHC isoform and at each MgATP condition examined (Figure 8). The distributions for the α-MyHC (Figure 8(a)) appear at shorter time periods compared to those of β-MyHC (Figure 8(b)). The distributions in both * * * * * * The two models agreed qualitatively in these respects, but did not agree quantitatively as estimates using the gamma distribution (GD) model were ∼40% higher than those using the single exponential distribution model (SED). (b) Values for t att using the two models correlated strongly. Despite the relative differences in the estimates of t att , either model produces similar results for comparing the effects of isoform and MgATP. * P < 0.05 by t-test between MyHC isoforms at same MgATP.
isoforms shift to longer time periods as MgATP is reduced. According to results from optical trapping [2,4], in both isoforms inorganic phosphate P i is released either prior to force production or very quickly after force production.  d) and 8(f)), the mean t att in both isoforms would have been prolonged due to longer rigor states. Particularly in the cases of lower MgATP, the value for t att would be made up of the sum of multiple times periods, as indicated by the increased shape parameter, and the distribution of t att is better represented by the gamma distribution.

Discussion
It is generally understood that force recorded in response to a length perturbation reflects (i) the number of myosin crossbridges bound to actin in a force-generating postpower stroke confirmation and (ii) a viscoelastic response that is the consequence of the length perturbation having been applied to temporarily formed crossbridges. The viscoelastic response, however, is often overlooked as a significant component of the dynamic force response. By frequency domain analysis, however, it is clear that the viscoelastic response is substantial in magnitude, and the frequency characteristics of the viscoelastic response emerge as a consequence of crossbridge kinetics [14,16,19]. The macroscopic measurements of viscoelastic mechanics in striated muscle combined with the modeling results presented in this paper can be used to estimate the microscopic temporal parameters reflecting myosin crossbridge kinetics, specifically the distribution of myosin attachment times (t att ) in a muscle fiber. Our current modeling culminated in a specific mathematic representation of the mechanical transfer function that would arise from intermittently attached myosin crossbridges, (9). This representation was also provided previously [19], but without significant discussion regarding the normalized survival function, that is, the PDF τ (t) term, in the integrand. We have provided here an explanation of the survival function, which underlies PDF τ (t) and its relationship with the distribution of t att , namely, PDF tatt (t) ((10a) and (10b)). Using the survival function and its relationship to the distribution of t att was a very important step in demonstrating how the macroscopic viscoelastic mechanics emerge from molecular phenomena. The gamma distribution model provided a more generalized representation of the many possible distributions of t att and in so doing provided an opportunity to discern the more subtle consequences of the t att distribution. The shape parameter in particular permits this more generalized description of the t att distribution. The effect of the shape The parameter values estimated for the α-MyHC demonstrate an increasing shape parameter and prolonged t att as MgATP was reduced. The resulting PDF tatt (t) looks more Gaussian in shape rather than exponential as MgATP was reduced. This would be consistent with a longer rigor state prior to MgATP binding and detachment of the myosin crossbridge. (b) The distributions for β-MyHC represent longer t att compared to α-MyHC. Again, a longer rigor state would arise with lower MgATP. (c, d) In the α-MyHC at 5 mM MgATP, the release of P i , release of ADP, and binding of ATP occur in relatively quick succession. The time periods that represent the P i and ADP states, however, would not be affected by MgATP availability. With the reduction in MgATP concentration, the rigor state would be prolonged as detected in the viscoelastic mechanics. (e, f). In the β-MyHC, the ADP state is longer than that in the α-MyHC and the rigor state is more dramatically prolonged with decreasing MgATP compared to α-MyHC.
parameter, which was always estimated to be near or greater than 2, indicates that the distributions of t att that underlie that recorded elastic and viscous moduli in both isoforms and under various MgATP conditions were distinct from the single exponential distribution and with lower MgATP begin to resemble Gaussian distributions. The increasing shape parameter with decreasing MgATP is consistent with the prolongation of the rigor state of the myosin crossbridge. We interpret these results for p and β to indicate that the crossbridge t att is best represented by two or more intermediate time periods, but we caution that a p approaching a value of 10 may simply reflect a longer t att and not additional biochemical states. The production of additional models stipulating only 2-or 3-specific states may be warranted. Such a model may be useful in detecting which biochemical state of the crossbridge cycle is affected by an intervention such as that due to drug or posttranslational modification. We are aware, however, that increasing the number of fitted parameters may also lead to greater covariance among the parameters and less meaningful results. We believe that our results demonstrate a benefit and utility of the gamma distribution model, which allows for near-zero occurrences of very short values for t att , as would be expected for some cases such as low MgATP concentrations and lower temperatures. The single exponential model, on the other hand, presumes that the most frequent values for t att occur at the shortest times. Our results, particularly at lower MgATP, suggest that the single exponential model may not be sufficient to fully describe the distributions of t att and the viscoelastic mechanical consequences as occuring in skinned muscle strips. It is notable, however, that while the gamma distribution model of t att provides a generalized model of the mechanical conseq uences of the intermittent attachment of myosin to actin, we found that the single exponential model is sufficient for purposes of comparing estimated values of mean t att based on the elastic and viscous moduli.
The gamma distribution model of myosin t att leads to characteristics of elastic and viscous moduli that might not be obvious or intuitively predictable. Specifically, we did not expect the fluctuations and asymmetry in the elastic and viscous moduli that emerged from a gamma distributed t att with p > 1. As the shape parameter increases and the t att distribution becomes more Gaussian in shape, the fluctuations and asymmetries observed in the elastic and viscous moduli in the frequency domain became more pronounced. Only through completion of the analytical derivation and computer modeling were we able to demonstrate these particular mechanical consequences.

Limitations.
Myosin kinetics are known to demonstrate dependencies on load or strain [10,12], which were not considered in the present work. The modeling in the present work focuses on explaining only the viscoelastic mechanics at higher frequencies and specifically ignores any strain dependencies thus providing a limited model of the macroscopic mechanical consequences of myosin kinetics. A more complete model of viscoelastic mechanics would need to predict and explain the dips in the elastic and viscous moduli like those observed at the lower frequencies. The negative viscous modulus, which corresponds to mechanical energy production by the muscle preparation, may well be one particular outcome of a strain dependency on myosin kinetics. For example, the myosin time attached may be shortened during muscle lengthening at these lower frequencies, perhaps in a P i -dependent manner [31][32][33], and the number of force-generating crossbridges attached during shortening versus lengthening would therefore contribute to an observation of mechanical work production.

Conclusion.
The importance of experimental observations of the single myosin molecule using optical techniques cannot be overstated. Without these methods, we would not likely be able to discern the nature of muscle force production at the molecular level. Discerning the performance of myosin within the context of myofilament lattice structure like that found in vivo is also important. Mathematical and computer models of muscle performance, such as that presented in this work, are valuable in their providing the opportunity to detect a measure of myosin performance, that is, the distribution of t att , within the context of in vivo sarcomeric structure that is not currently technologically possible by other means.

Appendix
In this appendix we provide a more detailed derivation of (15), which is the solution to (9) when PDF tatt (t) is represented by a gamma distribution. Therefore, PDF τ (t) in the integrand of (9) is represented by a normalized upper incomplete gamma function shown in (14): We choose to take the terms not containing t out of the integral: For purposes of solving the definite integral above, we scale every occurrence of t within the integral by β as follows: (A. 3) The definite integral has been solved previously [34]: Upon canceling factors we get the following: We choose to write the bracketed term in the following form: (A.6) Equation (A.6) represents the predicted complex modulus that would arise from an ensemble of myosin crossbridges whose time attached is described by a gamma distribution with shape parameter p and scale parameter β.
When p = 1 we get the special case for the gamma distribution being equivalent to a single exponential distribution. The solution provided in (A.6) for p = 1 then reduces to the following. The form of (A.7) is equivalent to that provided previously as a model for the C-process of sinusoidal analysis [19]. The mean time attached is represented here as β.