Molecular Mechanical Differences between Isoforms of Contractile Actin in the Presence of Isoforms of Smooth Muscle Tropomyosin

The proteins involved in smooth muscle's molecular contractile mechanism – the anti-parallel motion of actin and myosin filaments driven by myosin heads interacting with actin – are found as different isoforms. While their expression levels are altered in disease states, their relevance to the mechanical interaction of myosin with actin is not sufficiently understood. Here, we analyzed in vitro actin filament propulsion by smooth muscle myosin for -actin (A), -actin-tropomyosin- (A-Tm), -actin-tropomyosin- (A-Tm), -actin (A), -actin-tropomyosin- (A-Tm), and -actin-tropomoysin- (A-Tm). Actin sliding analysis with our specifically developed video analysis software followed by statistical assessment (Bootstrapped Principal Component Analysis) indicated that the in vitro motility of A, A, and A-Tm is not distinguishable. Compared to these three ‘baseline conditions’, statistically significant differences () were: A-Tm – actin sliding velocity increased 1.12-fold, A-Tm – motile fraction decreased to 0.96-fold, stop time elevated 1.6-fold, A-Tm – run time elevated 1.7-fold. We constructed a mathematical model, simulated actin sliding data, and adjusted the kinetic parameters so as to mimic the experimentally observed differences: A-Tm – myosin binding to actin, the main, and the secondary myosin power stroke are accelerated, A-Tm – mechanical coupling between myosins is stronger, A-Tm – the secondary power stroke is decelerated and mechanical coupling between myosins is weaker. In summary, our results explain the different regulatory effects that specific combinations of actin and smooth muscle tropomyosin have on smooth muscle actin-myosin interaction kinetics.


Video Analysis
Preprocessing of raw motility videos: Starting from raw motility videos, n consecutive frames are merged (brightness average of grey scale images). This lowers the number of frames to analyze, improves the signal to noise ratio by reducing brightness fluctuations, and reduces the influence of Brownian motion type noise (Fig. S1 A-C). In this study, every 10 frames were averaged into one, giving an effective time resolution of ∆t = 0.333 s. In each averaged frame, the brightness values are shifted by a scalar value so as to give the same brightness median across all frames to counter-act photobleaching and other lighting fluctuations. The brightness range of each video is rescaled to [0,1], so that the brightest pixel has the value 1, the darkest pixel 0.
Extraction of image objects: Each video frame is thresholded into a binary image, whose connected components are extracted (four neighbour connectivity). Objects of the size of one pixel are removed from further analysis right away.
Filament tracking: Filaments are tracked from one frame to the next based on a dissimilarity matrixM. For two consecutive frames (n − 1 and n), the dissimilarity matrix is an N n−1 × N n matrix, where N n is the number of filaments in frame n. The matrix elements are that is, the dissimilarity matrixM = (M k,l ) is determined from the area A k n−1 (number of pixels times area per pixel) and weighted centroidĉ k n−1 (center of mass of grey scale pixel values) of the k-th filament in the old (n−1) frame, and A l n andĉ k n of the l-th filament in the new (n) frame. Allĉ are two-dimensional vectors, and potentiation by 2 indicates the scalar product formed by the dot product of a vector with itself (i.e. the sum of squared elements).
The properties of a filament in a specific single frame n are stored under a Metaindex m. Given that in each frame n a Metaindex m uniquely references a single filament, the Metaindices can be used for filament tracking -a filament holding the same Metaindex m in two frames is understood as the same filament tracked over several frames.
The task of tracking, formulated in the above way, is to useM to assign Metaindices in a new frame n based on those from an old frame n − 1. First, the row minimum is found for all rows ofM. All values except these row minima are set to zero, e.g.
In this matrix, the three different possibilities by which a Metaindex can be assigned to a filament in the next frame are shown: (1) In the first and the second rows, the minimum is at the first position. This indicates that two Metaindices from the current frame point to the same filament in the next frame. In this case the filament in the next frame (in our example with index 1) is assigned a new Metaindex. (2) In the third row, the minimum is at the second position, so the Metaindex of the filament with index 3 in the current frame is forwarded to the filament with index 2 in the next frame. (3) The filament with index 3 in the next frame is not assigned a Metaindex from Metaindex forwarding, therefore it is assigned a new Metaindex. Whenever a filament in the next frame is not assigned a Metaindex, a new Metaindex is introduced. It is simply +1 greater than the greatest Metaindex existent at that point, which ensures uniqueness of Metaindices. Two more checks are executed during Metaindex forwarding: A length change check and a separation check. For the length change we compare the size between two filaments that are about to be assigned the same Metaindex in two consecutive frames. We check if 1) the ratio of the lengths of the n and n + 1 filament is greater than a defined parameter C rel or smaller than 1/C rel , or 2) the lengths of the n− 1 and the n filament differ by more than a defined parameter C abs . If so, we conclude that an unreasonably big length change occurred and assign a new Metaindex for the n filament. A separation check ensures that two filaments which were chosen to be assigned the same Metaindex in consecutive frames will actually form a connected trace. If they share at least one common pixel, the trace is connected, and the check is passed. If the filaments in consecutive frames do not share any pixels, they are assumed to be different filaments and a new Metaindex is introduced for the filament in the second frame. Similar approaches to following filaments were used in earlier works [1][2][3][4] Frame-to-frame velocity: The frame-to-frame velocity (V f 2f ) between two consecutive frames is the centroid displacement divided by the time resolution ∆t. V f 2f is used where time courses of filament motion are desired ( Fig. 1 B). However, for longer filaments an increasing underestimation of sliding velocity occurs due to the centroid cutting corners when the filament turns [5].
Filament length and trace-based sliding velocity: Filaments and traces, both of elongated shape, are transformed into rectangles of same area (A) and perimeter (P ). The object length then is the longer x + of the two rectangle edges [6] x The trace-based velocity (V ) is calculated from the distance that a filament travelled along its trace during the time it was observed and tracked for (T ). The trace is constructed by merging the filament's binary images from all frames ( Fig. 1 B). The filament length (L) is determined by averaging over the filament length in all frames that the filament was observed in. L is subtracted from the trace length, resulting in the distance travelled by the filaments tip along its sliding path ( Fig. 1 B). This trace length is divided by T to give V . Unlike the centroid method, V is not affected by systematic errors resulting from curved motion.
Features of actin motility: The mean sliding velocity is calculated as the mean of trace velocities V . Trace velocities are used instead of the centroid-based frame-to-frame sliding velocities (V f 2f ) to prevent velocity under-estimation for long filaments ( Fig. 1 B) [5]. The motile fraction (f mot ) is the fraction of V f 2f greater than a threshold velocity V t (Fig. 1 C). Based on another threshold velocity . Periods spent in the non-motile state without transition into the motile state are called stop times t stop , conversely periods spent in the motile state without transition into the non-motile state are called run times t run (Fig. 1 C). In this study, V t = 0.125 µm/s and V t t = 0.225 µm/s were used.

Assessment of video analysis
Influence of random displacement: Filament centroids get displaced due to directed filament motion and random fluctuations, which were approximated as a Brownian-motion type influence on V f 2f . The mean velocity of directed sliding (ν) is independent of the time resolution ∆t. The velocity from Brownian motion is proportional to (∆t) −1/2 ( Fig. S1 A, B). This implies that excessively lowering ∆t will obscure V f 2f distribution features which are visible at sufficiently high ∆t (Fig. S1 C). For sufficient ∆t, a motile and a non-motile population in the V f 2f can be distinguished [7], indicating that a sufficiently large ∆t was chosen to discern these features of in vitro motility. Accuracy of L and V values: Computer-generated mock motility videos were used to assess the analysis's accuracy. V and L were faithfully detected above the diffraction limit. Below the diffraction limit, indicated by complex L or V (Eq. 3), V and L become quantitatively unreliable, but filaments are still qualitatively ordered by L (Fig. S2 A). At a too coarse time resolution, L is overestimated as filament motion over several consecutive frames "stretches" filament images similar to a motion blur (  [8]. However, less of the filaments successfully pass the rejection criteria (Fig. S2 C).
L below the diffraction limit: In the rectangular transformation to determine filament and trace lengths (Eq. 3) complex solutions with an imaginary component = 0 arise for P/4 > √ A. This is the case where filaments that are shorter than the diffraction limit of ≈ 0.2µm appear in the videos as approximately circular objects ( Fig. S2 A-C). Below the diffraction limit filament and trace lengths are not quantitatively precise. Still, longer filaments can be expected to hold a greater number of fluorophores, increasing the brightness and in turn the optical area of the filament, resulting in a greater image object. In consequence, the real parts P/4 could be used to qualitatively sort affected filaments according to L, even though quantitatively accurate L could not be inferred (Fig. S2 A-C).
Relevance of filament width and digitization threshold: Based on analysis results obtained with different Black-White thresholds (BW ), V and L averages depend on filament width (W ) (Fig. S3). However, in a Black-White threshold range between 0.2 and 0.35, V , W , and the number of analyzed filaments seem mostly unaffected. For a low BW < 0.15, the computation time increases sharply, indicating that optical noise and brightness fluctuations are analyzed as objects (Fig. S3, inset). L depends on BW across all BW values. This can be understood when taking into consideration that a lower BW leads to an increase in the visual area covered by objects, giving a "blobbier" appearance. In consequence, also the head and tail tips of filaments reach out further, effectively increasing L. For BW > 0.55 the number of valid filaments drops sharply, indicating that an excessively high BW effectively removes most filaments.

Machine Learning Quality Control
Filament traces can be rendered unfit for analysis by filament crossings, filaments "swirling" around attachment points on the motility surface, lateral displacement [9], or other features that lead to inaccurate V and L estimation. Initially, such traces were removed manually, leading to acceptance (1) or rejection (0) of specific traces. Decision tree models were then trained on such manually generated data sets, to replace time-consuming manual inspection of filament traces. To make trace images more uniform and thereby easier to categorize, trace images were rotated so that all images' main axes pointed in vertical direction. The binary image data; length, width, solidity, perimeter, and area of the trace; and the number and position of corners found using MatLab's corner function were extracted as features. These features were used to train the decision tree models on manually scored data. Using > 100 decision trees did not lead to a further decrease in the prediction error rate, indicating ≈ 100 as the number of decision trees necessary to achieve the maximally possible classification performance (Fig. S4 A). Models used in this work were based on 150 decision trees, ensuring a sufficient number of decision trees. Receiver Operating Curves for the classification of data from one day based on training on data from another day show that a good margin for successful classification exists (Fig. S4 B). Based on cost optimization for different data set combinations, we chose an acceptance quorum of 0.7 trees accepting a filament trace (Fig. S4 C). False positives (traces that should have been rejected) were visually inspected and showed little or no features deviating from an undisturbed, elongated shape (data not shown).

Statistical Analysis
The raw data are motility assay videos recorded for the different conditions (Tab. 1), in the following the conditions are indexed by C. In each condition, a number (V C ) of different videos were recorded, in the following the videos are indexed by v = 1, . . . , V C . Let v be a sequence of length V C , whose elements each refer to one video of C. In a general case, a video can occur in v once, more than once, or also not at all. The meaning of v is to describe a data set made up of the different videos referenced in v, e.g. v = (1, 2, . . . , V C ) represents the original data set in which each video of condition C is referenced exactly once. Resolution by L is executed using a sliding window approach: an overall L range (L ∈ [L min , L max ]) is set, then the center of a window of constant width is moved along this range at equally sized steps. At each step, all filaments from the videos contained in v (multiply referenced videos are included multiply) within the windows current L range are gathered, the specific motility feature is calculated for all filaments, and the arithmetic mean of these feature values is stored as the according element of a result vector (R C v ). In this study, 50 length windows were applied for each of the four features (Fig. 2), resulting in an R C v with a length of 200 elements. Within this formalism, R C v=(1,2,...,V C ) refers to the original data set with all videos contained exactly once; this can be likened to a "mean" in a scalar measurement. However, when the videos referenced by v are assigned by random draws with replacement, effectively bootstrap resampled R C v are created, which serve as empirical distributions around the "mean". Statistical assessment in Principal Component space: PCA was applied to reduce feature dimension and remove correlations between different elements in the R C (v omitted from notation) values to prevent inflation of statistical significance. First, the bootstrap R C from all C were pooled to execute the PCA on all bootstrap data at once. The elements of the bootstrapped R C were transformed into z-scores (centered on 0 by subtraction of the mean and normalized by standard deviation in each element) and the first 3 PCs were considered in the statistical assessment (Fig. 3 A, B).
After transformation into PC space, the overlap of 95% intervals of the empirical distributions of the bootstrapped R C was assessed to detect if statistically significant differences with p < 0.05 are present. To compare two conditions C = I and C = II, the original (not bootstrapped) R I and R II PC coordinates were determined using the PCA transformation matrix determined based upon the bootstrapped R C , giving r I and r II , respectively. A vector d = r II − r I connecting both C = I, II was determined and the bootstrapped R I and R II were projected onto that vector using a simple scalar (dot or inner) product. Then, empirical 2.5% and 97.5% percentiles for each C were constructed from these one-dimensional values (Ward's bootstrap confidence intervals). Where no overlap between these percentiles existed, statistical significance was inferred (p < 0.05) (Fig. 3 C, D).
While PCA and bootstrapping are standard methods, their combined application is seen more rarely. In other studies, for example, bootstrapping is used to estimate variance in the PCA procedure itself [10,11], or a data set of several measurements is bootstrapped into mock data sets for each of which an independent PCA is executed to determine statistics of the PCA scores [12,13].
Hierarchical clustering: The bootstrapped R C in PC space were pooled across all C, and then subjected to agglomerative hierarchical clustering to rediscover the conditions (Fig. 3 E-H). The Ward method and an Euclidean metric were used.
Length-averaged fold change analysis: αA was chosen as a reference condition (C = 0), by which the other conditions were normalized. All bootstrapped R C are divided element-wise by the values of the reference condition's R 0 original (not bootstrapped) feature vector s C values. From the normalized bootstrapped R C , the V , f mot , t stop , and t run are then averaged over L (arithmetic mean), resulting in scalar fold change values.
To estimate the rate of type I errors (false positive rate), the rate of rejection of the null hypothesis in absence of an actual difference was determined. We executed length-averaged fold change analyses as described above, with the difference of comparing the baseline condition (αA) to itself. Due to the bootstrapping employed in this procedure, two executions of the procedure correspond to a random reshuffling of the same data set. This preserves the same original data set, thus creating two data sets that have no difference between them but do account for the fluctuations being present in the recorded data. The rejection of the null hypothesis, i.e. the detection of statistically significant differences, applied when the confidence intervals did not overlap. This would be the case if where CI i,j are the limits of the confidence intervals. i ∈ {1, 2} refers to the first and second resample, respectively. j ∈ {1, 2} refers to the smaller and greater limit of the confidence interval, respectively. In consequence, a significant difference would be detected if Empirical distributions of ∆ CI indicate that, given our specific assessment method, number of measurements, and fluctuations in our data, the detection of statistically significant differences in the absence of actual differences is highly unlikely (Fig. S5).

Mathematical Model and Simulation
We use a mathematical model developed in one of our earlier studies [14]. Here, we give only the essential assumptions, expressions, and descriptions. We simulate N myosin binding sites (N = L/0.0355 µm [15,16] This means that only the main power stroke of 4 nm length affects the equilibrium position, while the secondary power stroke of 2 nm is followed by an immediate detachment (ATP-saturated buffer) [17]. The kinetics of the individual binding sites consist of a unidirectional reaction cycle (0 → 1 → 2 → 0), with transition rates k a (attachment, 0 → 1), k p (main power stroke, 1 → 2), and k d (detachment, 2 → 0) [17].
Due to actin acting as a rigid backbone, the x m result from an interaction between all myosins currently bound to actin (σ m = 1). The length travelled by the actin filament (a, in units of µm) also results from this mechanical interaction. For a binding site going from c m = 0 to c m = 1, an unstrained position is assigned, which is a random variable drawn from a normal distribution around the mean a (attachment without any strain) with a standard deviation of w. The force exerted by myosin on an individual binding site m is K is a spring constant, we assume linear elasticity. After each kinetic transition, we assume that an instantaneous force equilibrium between all attached myosins is attained, a is left unchanged if n att = 0. Transition rate calculation: The current chemical state c m as well as the work that has to be exerted for the next kinetic transition (∆W m ) determine the rate at which a binding site m will undergo the next transition, r m : 2c f /K is a coefficient quantifying how strongly mechanical work exerted during a mechanical transition affects the rate of this mechanical transition. When a, x 0 m and c m are known, the overall mechanical work stored in the actin-myosin system can be calculated: m denotes the general summation over all myosin binding sites, irrespective of the current value of m. We now assume that a is continuously changed to the system's mechanical equilibrium position throughout the execution of a mechanical step. It follows that ∆W m = W af ter m − W bef ore m (13) can be calculated from the work before (W bef ore m ) and after (W af ter m ) a mechanical step of size d occurring at binding site m, The simulation is initialized with a = 0, x m = 0, and randomly assigned c m . Filament motion is simulated by iteration of the following steps.
1. Instantaneous force equilibration by adjustment of a.
2. Current transition rates r m are calculated based on the system's current state.
3. The next transition time and the binding site undergoing a transition is stochastically drawn based on a Gillespie scheme [18]. The according binding site is advanced by one step in the reaction cycle.
This sequence is repeated until either a maximal number of kinetic transitions per binding site (n steps ) or a maximal simulation time (t max ) is reached. The iteration results in a sequence of travelled lengths a, and elapsed times t. In the actual video analysis from the motility assay, these exact time courses are not accessible, instead only the difference in travelled length at time intervals ∆t (time between consecutive frames in our videos) can be assessed. Accordingly, V f 2f were extracted from the stochastic simulation using an imposed time resolution of ∆t. Actin sliding was simulated for N = 5, 6, . . . , 90. At each N , actin sliding was simulated till 100 s of simulation time or 100, 000 · N chemical transitions were reached. Features of actin motility were extracted for each N , followed by a running window analysis (Start: N = 10, end: N = 90, window width of 20). Considering a binding site distance of ≈ 0.035 µm [15,16], this approximates the L range we analyzed in our experiments.
Correction for non-specific binding to motility surface: When this model was simulated, above N ≈ 25, t run and t stop could not be measured due to V f 2f not falling below the threshold for determining t run and t stop . Occasional non-specific stopping events occur for longer filaments [14]. We introduce these stopping events by inserting Poisson-distributed V f 2f = 0 µm/s values into the V f 2f time trace: where ∆T is the time between two consecutive stop events, r a random number drawn from a uniform distribution r ∈ [0, 1]. T is the mean time between two consecutive stop events, which has been adjusted to limit t run as observed in our experiments (Fig. 2 D).
Parameter adjustment: The baseline model parameters (Fig. 7) were taken from the literature and then adjusted to fit the L resolved V and f mot (k a [19,20], k p [17], k d [17], an additional factor of 5.5 was used to account for higher temperature in our experiment [14]), or estimated from L resolved V and f mot , starting from values determined in our earlier study of skeletal muscle myosin [14] (c f , w). The altered parameters mimicking the αA-Tmαβ, γA-Tmβ, and αA-Tmβ conditions were estimated based on Fig. S6 and comparison with L resolved features and scalar fold changes in the conditions with statistically significant differences in the experiment ( Fig. 2 and 3, respectively).