Direct Measurements of Local Coupling between Myosin Molecules Are Consistent with a Model of Muscle Activation

Muscle contracts due to ATP-dependent interactions of myosin motors with thin filaments composed of the proteins actin, troponin, and tropomyosin. Contraction is initiated when calcium binds to troponin, which changes conformation and displaces tropomyosin, a filamentous protein that wraps around the actin filament, thereby exposing myosin binding sites on actin. Myosin motors interact with each other indirectly via tropomyosin, since myosin binding to actin locally displaces tropomyosin and thereby facilitates binding of nearby myosin. Defining and modeling this local coupling between myosin motors is an open problem in muscle modeling and, more broadly, a requirement to understanding the connection between muscle contraction at the molecular and macro scale. It is challenging to directly observe this coupling, and such measurements have only recently been made. Analysis of these data suggests that two myosin heads are required to activate the thin filament. This result contrasts with a theoretical model, which reproduces several indirect measurements of coupling between myosin, that assumes a single myosin head can activate the thin filament. To understand this apparent discrepancy, we incorporated the model into stochastic simulations of the experiments, which generated simulated data that were then analyzed identically to the experimental measurements. By varying a single parameter, good agreement between simulation and experiment was established. The conclusion that two myosin molecules are required to activate the thin filament arises from an assumption, made during data analysis, that the intensity of the fluorescent tags attached to myosin varies depending on experimental condition. We provide an alternative explanation that reconciles theory and experiment without assuming that the intensity of the fluorescent tags varies.


Parameter estimation
In order to model the measurements of Desai et al. [4], we had to estimate four parameters: σ N the standard deviation of the assumed Gaussian background noise; σ F the standard deviation of the assumed Gaussian temporal fluctuation in GFP intensity; e the fluorescent emission of a GFP; and ε which determines the degree of regulation. All parameters, except ε, were estimated prior to fitting the data with the model. We varied ε in order to optimize the fit of the model to the data.
We estimated the background noise, σ N , individually for each experimental condition. The values used in the model are listed in Table 1 of the main text. To estimate σ N from a kymograph, we first removed spatial fluctuations in intensity. To do so, we calculated the minimum fluorescence at every pixel along the actin filament. We then fit a fifth-order polynomial to this minimum fluorescence as a function of position along actin. The resulting curve defined zero fluorescence for every pixel along actin.
We then plotted a histogram of the fluorescence of every pixel in a kymograph. The data show a peak, which corresponds to the background noise, and an extended tail toward positive fluorescence, which corresponds to the binding of fluorescent myosin (GFP-S1). We estimated the standard deviation of the histogram by matching a Gaussian to the distribution to the left of the mean, thereby only fitting the background noise (see Fig. 1A).
To estimate σ F and e, we analyzed a histogram, which plots the frequency of spots of a given intensity, measured under conditions where we expected to see primarily single molecules binding (pCa 6, [ATP] = 0.1 µM, [Myo] = 1 nM). The maximum of this peak was around I 1 = e/f = 45 intensity units (Fig. 1B). Given that these data were collected at a frame rate of f = 10 Hz, this gives a value of e = f I 1 = 450s −1 .
Since the histogram is constructed by fitting the raw data with a Gaussian, we expect that the contribution of σ N to the signal is small. One might then expect the histogram to be Gaussian with standard deviation σ F , but the histogram is not exactly Gaussian. This is likely due to imperfections in the fitting algorithm used to construct the histogram. However, the histogram is approximately Gaussian near the peak, so we estimated σ F = 0.22 from matching a Gaussian to this peak. In support of the view that the non-Gaussian shape of the histogram is due to the fitting algorithm, and that the width of the peak can be reasonably used to estimate σ N , simulations with only single binding events generate histograms that are similar to observations (Fig. 1B). Hz. The spread in the histogram near the peak is due to temporal fluctuations in intensity, and is well-fit by a Gaussian with standard deviation σ F = 0.22 (blue), in units of scaled intensity. Deviations from the Gaussian are likely due to due to the algorithm used to determine the fluorescent intensity of a spot, since simulations with only single binding events generate histograms that are similar to observations (gray).
To estimate ε, we fit the data [4]. A series of preliminary simulations suggested that ε ≈ 0.06 fit the data at pCa 6 the best. To quantify this, we performed simulations of the experiments at variable [Myo] (pCa 6, ATP = 0.1µM, [Myo] = 1, 5, 10, 15nM, 1000 frames collected at 10Hz), with ε = 0.04, 0.06, and 0.08. We compared both the distributions of clusters ( Fig. 2A), and mean fluorescence per pixel (Fig. 2B). In all cases, ε = 0.06 was the best (Fig. 2C). We therefore used this value for pCa 6. Since there was only a single measurement at pCa 5 and pCa 7, we did not perform as detailed an analysis, but rather estimated ε = 0.01 and ε = 0.4, respectively, from preliminary simulations.

Average fluorescence
In the main text, the comparison between model and data is presented qualitatively in Figure 4A-C, where simulated and measured kymographs are displayed side-by-side. To obtain a more quantitative comparison, we determined the average fluorescence per pixel in each kymograph. In the main text, we present only measurements at variable [Myo] (Fig. 4D). All measurements are shown in Fig. 3. The agreement between simulation and measurement is good.
3 More details of constant vs. variable GFP emission 3

.1 Desai et al's analysis
The model of local coupling between myosin molecules [1,2,3] assumes that a single myosin can activate the thin filament; Desai et al [4] conclude that two myosin heads are required to activate the thin filament. Desai  Figure 3: Plots of average pixel fluorescence for kymographs measured under various conditions. The left plot, at variable myosin, can be found in the main text as Fig. 4D. The other two plots show variable ATP and variable calcium. The model (solid dots, SD error) reasonably agrees with the data from Desai et al. [4] (hollow dots). Note that the lowest ATP measurement (middle plot) was simulated at [ATP]=0.15 µM, as discussed in the main text. For each pixel, zero fluorescence is defined as the minimum value it achieves during the kymograph. In each plot, the apparent background fluorescence is indicated as a dashed line. Fluorescence is measured in scaled intensity, defined in the text, which is non-zero in the absence of a myosin, and a single myosin increases the signal by 1 unit. et al's [4] conclusion follows from their analysis of the histograms they measured (see Fig. 2 of the main text for the construction of these histograms, and Fig. 5 of the main text for the histograms themselves). We now summarize this analysis.
Each measured histogram contains the fluorescent intensity of every spot in a 500 or 1000 frame movie. Each fluorescent spot comes from a cluster of GFP-tagged myosin molecules (GFP-S1s). Given that 1. each GFP-S1 has an average fluorescent intensity I 1 2. a cluster of i GFP-S1s has a measured intensity of iI 1 ± σ S , where σ S represents signal noise 3. signal noise is Gaussian and its magnitude is independent of the number of GFP-S1s in the cluster then each histogram F (I) arises from the following sum The coefficients, a i , determine the relative frequency of clusters of i molecules in the histogram. Thus, for example, if a histogram consists only of clusters of single GFP-S1s, then a 1 = 1 and a 2 = a 3 = a 4 = · · · = 0. Desai et al [4] determined the coefficients, a i , by first measuring σ N for an isolated GFP and then fitting their histogram with the following equation allowing the a i and I i 's to vary in order to optimize the fit (Fig. 4A). Each time they performed this fit, they found that the optimal I i values occurred at regularly spaced intervals -consistent with Eq. 1, where I i = iI 1 . Thus, one might expect that the lowest I i is I 1 , corresponding to a cluster of one GFP-S1, the next 2I 1 , corresponding to a cluster of two GFP-S1s, and so on. If so, then if the I i 's are sorted from small to large, and plotted as a function of apparent cluster size (i.e. 1 for the first, 2 for the second, and so on), the resulting curve should be linear with slope I 1 and should extrapolate to 0 at a cluster size of 0. Although linear curves were always observed, frequently the line did not extrapolate to 0 (Fig. 4B, black curve). Instead, the curve extrapolated to 0 only upon assuming that the first cluster contained two GFP-S1s, the second three, and so on (Fig. 4B, red curve). Based on this observation, Desai et al [4] conclude that single GFP-S1s cannot bind to the thin filament; rather, two or more GFP-S1s are required to activate the thin filament.

Analysis of simulated data
Given that the model, which assumes a single GFP-S1 can activate the thin filament, successfully reproduces the measurements, one might guess that the measured histograms do not necessarily imply that two or more GFP-S1s are required to activate the thin filament. However, there are subtle differences between the measured and simulated histograms. It is therefore possible that the analysis of Desai et al [4] is sensitive to these subtle differences, so that the simulated data, when analyzed, will differ from the analysis of the measurements. We therefore performed Desai et al's [4] analysis on the simulated measurements with variable [Myo] (Fig. 4). We used a slightly different fitting algorithm than Desai et al [4]. Our algorithm is as follows: 1. Use linear interpolation to represent a given histogram with 1000 equally spaced points between 0 and 600.
2. Determine the intensity (I max ) at which the maximum frequency (F max ) occurs.
3. Use a non-linear optimization algorithm (Matlab's fminsearch function), starting from an initial guess of mean I max and amplitude F max , to fit the histogram with a single Gaussian, with standard deviation σ N , allowing its mean and amplitude to vary to minimize the least squares error between the Gaussian and the interpolated histogram 4. Determine the intensity (I E ) at which the measurement and fit have a maximum difference (F E ). A. In the first step of the analysis, a histogram is fit with Gaussians of fixed standard deviation. Here, a simulated histogram (pCa 6, ATP = 0.1µM, [Myo] = 15nM, 1000 frames collected at 10Hz) is fit with eight Gaussians, each of standard deviation σ = 15. B. In the second step of the analysis, the Gaussians are ordered, numbered sequentially and plotted. In all cases, a linear curve results. Sometimes, that linear curve does not pass through the origin (black), but when the first point is assigned to two molecules, the linear curve passes through the origin (red). Note: the starred point is considered an outlier and is not included in the linear fit. C. In the third step of the analysis, the amplitudes of the Gaussians that correspond to each number of binders are plotted. We wrote our own algorithm for this analysis, which differs slightly from Desai et al's [4] algorithm. We used our algorithm to analyze measured histograms (blue), and simulated histograms (gray, individual, black, average, yellow, range) collected at 10Hz, variable myosin, pCa 6, ATP = 0.1µM, 1000 frames. Analysis of the measurements and simulations with our algorithm yields similar results to those reported in Desai et al. [4] (blue squares). 5. Use a non-linear optimization algorithm (Matlab's fminsearch function), starting from an initial guess of the previous best-fit with one additional Gaussian, with standard deviation σ N , of mean I E and amplitude F E .
6. Repeat steps 4 and 5 until no improvement occurs, and/or two Gaussians with similar means are observed.
To validate this algorithm, we analyzed the measured histograms at variable [Myo], allowing σ N to vary in order to optimize the fit to Desai et al's [4] analysis of these histograms. The results, shown in Fig. 4C with Desai et al's [4] analysis as blue squares and our fit as blue lines, are in reasonable agreement. Further, the best-fit values of σ N = 9, 12, 17 and 15 at [Myo] = 1, 5, 10 and 15 nM, respectively, are in reasonable agreement with our estimates (see Parameter Estimation in Supplementary Material) of σ N = 9.5, 11.25, 13.5 and 16.2 at [Myo] = 1, 5, 10 and 15 nM, respectively.
When we used this algorithm to analyze the simulated data, the results are similar to those from the measurements (Fig. 4C). Further, many of the fits exhibited the same non-zero extrapolation reported in Desai et al [4] (e.g. Fig. 4B is from a simulation at [Myo] = 15 nM). While these results are the basis for the conclusion that two or more myosin molecules are required to activate the thin filament, these simulations assume that a single myosin can activate the thin filament.