A Python based automated tracking routine for myosin II filaments
AA Python based automated tracking routine formyosin II filaments
L.S. Mosby , , M. Polin , , D.V. K¨oster Centre for Mechanochemical Cell Biology, University of Warwick, Coventry CV47AL, UK Physics Department, University of Warwick, Coventry CV4 7AL, UKE-mail: [email protected], [email protected]
Abstract.
The study of motor protein dynamics within cytoskeletal networks isof high interest to physicists and biologists to understand how the dynamics andproperties of individual motors lead to cooperative effects and control of overallnetwork behavior. Here, we report a method to detect and track muscular myosinII filaments within an actin network tethered to supported lipid bilayers. Based on thecharacteristic shape of myosin II filaments, this automated tracking routine allowed usto follow the position and orientation of myosin II filaments over time, and to reliablyclassify their dynamics into segments of diffusive and processive motion based on theanalysis of displacements and angular changes between time steps. This automated,high throughput method will allow scientists to efficiently analyze motor dynamics indifferent conditions, and will grant access to more detailed information than providedby common tracking methods, without any need for time consuming manual trackingor generation of kymographs. single particle tracking, automated detection, myosin
Submitted to:
J. Phys. D: Appl. Phys. a r X i v : . [ q - b i o . S C ] J a n fficient tracking of myosin II filament position and orientation over time
1. Introduction
Molecular motors are important for many cellular processes such as cell cortexdynamics, cell migration, and the intracellular transport of vesicles. Purification ofvarious molecular motors from tissue samples or after recombinant expression, andthe analysis of their biochemical and bio-physical properties in reconstituted systems,was instrumental in furthering our understanding of the motors and how they areregulated in live cells. Whereas most molecular motors operate individually or in dimers,myosin II motors usually form bundles of about 10 (non-muscular myosin II) up to 250(muscular myosin II) proteins at physiological salt conditions, giving rise to cooperativeeffects between the myosin head domains that govern the effective binding dynamics,speed and processivity of these myosin II protein ensembles. After these cooperativeeffects were studied theoretically, recent advances in microscopy and the design ofreconstituted acto-myosin networks have allowed the study of the dynamics of myosinII filaments experimentally [1, 2]. Traditionally, motor dynamics were studied usingkymograph analysis or common single particle tracking (SPT) routines that identifypoint-like structures. These methods have drawbacks, for example kymographs onlytrack dynamics along a given path and thus cannot be used to study two-dimensionaldiffusive motion. Similarly, common SPT routines cannot take into account the myosinfilament orientation, which is a valuable parameter that can be used to characterisetheir dynamics. Here, we applied an image analysis routine from Astrophysics thatwas designed to identify galaxies [3–5] to detect myosin II filaments, as both types ofobject are characterised by their elongated, elliptical shapes with limited signal to noiseratios in large, heterogeneous samples. After particle detection, tracks are generated bycomparing the positions, orientations and detected areas of particles between subsequentframes. This Python based routine worked robustly in the detection of myosin IIfilaments in a large, crowded sample set and made it possible to distinguish betweenthe diffusive and processive motion of the filaments.
2. Methods
The experimental data used here to demonstrate the detection and tracking ofmyosin II filaments was obtained by interferometric scattering microscopy (iSCAT),as described previously [2, 6]. Actin filaments and skeletal myosin II filaments werepurified from chicken breast muscle following established protocols. Glass coverslips( . before being used for the preparationof experimental chambers. After formation of supported lipid bilayers (containing98% DOPC and 2% DGS-NTA(Ni ) lipids (Avanti Polar Lipids Inc., US)) andaddition of our actin-membrane linker protein HKE (decahistidine-ezrin actin binding fficient tracking of myosin II filament position and orientation over time mM KCl, mM M gCl , mM EGT A, mM HEP ES, pH .
2) to allow the formationof membrane tethered acto-myosin networks [2, 7]. Images were recorded on an iSCATmicroscope setup equipped with a cMOS camera and processed to reflect interferometriccontrast values in 32bit with a pixel size of 0 .
034 x 0 . µ m [2]. We used the Python programming language as a platform for the myosin II filamentdetection based on the Python library for Source Extraction and Photometry (SEP)[3–5], which generates the position, spatial extent, and orientation of ellipsoidal particlesfor each frame of a time lapse image sequence.To distinguish individual myosin II filaments, we first applied a signal and areathreshold to isolate probable myosin II filament detection from background, followed bythe application of a refined area threshold to distinguish individual myosin II filamentsfrom aggregates and to reduce the rate of false-positive detection due to poor signalto noise ratios. For the first step, a detection was defined as a region of >
40 pixels,each with an intensity of at least 1 . σ g,rms above the local background intensity, where σ g,rms is the global root-mean-square error of the spatially varying background of theimage. These minimum area and intensity cutoffs allowed confident detection of myosinII filaments without falsely detecting background fluctuations.The SEP Python package computes a set of ellipse parameters (a: semi-major axis,b: semi-minor axis, θ : orientation on a − π/ ≤ θ ≤ π/ A = πab (figure 1(a,b)). For clarity the ellipse parameters were scaled to x = 6 x detected ( x = a or b ) to allow for better representation of the detected object byeye [5]. It should also be noted that the observed signal is a result of the convolutionbetween the real signal and the point-spread function of the microscope. For thesereasons the areas calculated using this method have only been used for diagnostics.In order to prevent the false-positive detection of background intensity fluctuations(figure 1(c)) a minimum area cutoff was implemented. After considering the probabilitydistribution function of the detected myosin II filament areas for a sample data setwith a high signal to noise ratio, the minimum area cutoff was defined as equal tothe smallest detected area A min = 0 . µ m . An average area for a single myosinII filament detection was also calculated as A av = 1 . A mode, where A mode wasthe maximum value of an exponentially-modified Gaussian distribution fitted to thearea probability distribution function of a data set recorded at high ATP concentrationdisplaying minimal myosin II filament aggregate formation (figure 1(e)). At lowerATP concentrations the sample data sets contained many myosin II filament aggregates(figure 1(d)), which resulted in a higher than expected rate of detection of large areavalues contributing to the exponentially decaying tail of the area distribution shown in fficient tracking of myosin II filament position and orientation over time A min, equalto the minimum area detection in the high signal-to-noise ratio data set with the lowestATP concentration. (b) Log-log plots of the distributions showing their decay profilescompared to lines of constant gradient in log-log space. (c) Example of anomalouslysmall area detection in the data set with the highest ATP concentration (red) due to poorsignal to noise ratio and a disperse region of pixels with intensity greater than 1 . σ g,rmsabove the local background intensity. (d) Examples of myosin II filament aggregates inthe data set with the lowest ATP concentration (red). (e) The corrected distributionafter the removal of anomalously small area detection from the data set with the highestATP concentration. An exponentially-modified Gaussian distribution was fitted to thedata (orange), and the red vertical line shows the average area A av = 1 . A mode.figure 1(a,b). In order to observe the dynamics of single myosin II filaments, tracksincluding larger aggregates must be removed before analysis. fficient tracking of myosin II filament position and orientation over time Tracks of myosin II filament motion were generated by considering the spatialdisplacement and the change in area of the detected particle between subsequentframes. First, it was checked whether the change in position of a detection between twoconsecutive frames was suitably small. An estimation of the maximum bound velocityof a moysin II filament of v max ∼ . µ m s − suggests filaments can move a maximumdistance of | ∆ x | max ∼ | ∆ x | = (cid:112) (¯ x − ¯ x ) + (¯ y − ¯ y ) ≤ x, ¯ y are the co-ordinatesof the centroid of the detected particle, and the subscripts indicate initial (1) and final(2) positions), to allow for small deviations in motion and shape.Second, a detection at time t was also only added to the track of a detection attime t if its area A ( t ) satisfied A ( t ) − A av ≤ A ( t ) ≤ A ( t ) + A av, which allowedfor the transient overlap of two detections. This can lead to tracks being cut intomultiple shorter tracks due to overlap events, which would skew dwell time and totaldisplacement distributions towards lower values. This effect is particularly noticeable indata sets with lots of aggregation or clustering of particles. Based on the intensity andarea thresholds, it can be assumed that the detected particles must have been boundfor the majority of the exposure time (here 200 ms). However, for the remaining areaanalysis a detection only associated with tracks of at least two frames duration areincluded, to further prevent noise affecting the result.In order to derive an appropriate maximum area threshold that minimises thedetection of myosin II filament aggregates, the effect of varying the maximum area (asa multiple of the average area, A av) on average system variables was analysed (figure2(a,b)). The maximum area was set to A max = 4 A av = 3 . µ m as it maximisedthe measured average dwell time, while also being larger than the value that maximisedthe measured average distance travelled. It should be noted that all tracks that containa point with an area A > A max have been removed in their entirety from the datato be analysed in order to prevent artificially skewing the dwell time and displacementdistributions. In order to measure dwell times and total displacements as accurately aspossible, only tracks that ended (by detachment or by leaving the field of view) by theend of the video were included in the analysis. The effect of myosin II filaments leavingthe field of view on the measured dwell time was found to be negligible, as less than4% of tracks ended within 6 pxl (the maximum allowed filament displacement betweentimesteps) of the edge of the video. The final probability distribution function describingthe detected areas after applying the above corrections is shown in figure 2(c,d).When two myosin II filaments transiently overlapped, the aggregate was alwaysappended to the track with the current longest dwell time in order to probe longertrack dynamics. The ratio of tracks that ended within l = √ A av of any point onanother track (potentially due to aggregation) was independent of the dwell time ofthe tracks (results not shown). Similarly, the dwell time distributions and associated fficient tracking of myosin II filament position and orientation over time A max (in multiples of A av) on the average dwell time of the myosin II filaments in the iSCAT videos forthe highest and lowest ATP concentration data sets. (b) The effect of varying themaximum area cutoff A max on the average straight-line distance the myosin II filamentstravelled between their initial and final positions in the iSCAT videos. (c) The correctedprobability distribution functions of detection areas after the application of the derivedminimum and maximum area cutoffs. Red vertical line shows the maximum area cutoff.(d) Log-log plots of the corrected area distributions showing their more favourable decayprofiles.characteristic timescales, including the average dwell time, of the tracked myosin IIfilaments at varying ATP concentrations did not change significantly after the removalof the tracks that ended in an overlap event (results not shown). These results suggestthat the overlap events do not introduce a characteristic timescale into the dwell timeanalysis, and hence can be ignored.Following particle detection and tracking, analysis is carried out as follows: fficient tracking of myosin II filament position and orientation over time Of particular interest in this work is the splitting of tracks into regions of processive,directed motion and diffusive, non-directed motion. The processive, directed motion ofthe myosin II filaments executed when the filaments are fully bound to the underlyingactin network has a much longer persistence length than the diffusive motion exhibitedby the filaments when they are only partially bound to the surface. Tracks were definedas being diffusive or processive using the formalism from the work by Jeanneret et al. [8];firstly tracks that contain large, directional displacements that cannot be a result ofBrownian motion were identified, and then the exact frames corresponding to thesedisplacements were isolated by considering the correlation in the displacements betweenadjacent timesteps.A purely diffusive, spherical particle with a time-dependent position x ( t ) anddiffusivity D in two-dimensions will always have an average displacement of zero anda mean-squared displacement equal to (cid:104) ∆ x ( t ) (cid:105) D = 4 Dt . Assuming that the myosinII filaments being tracked exhibit non-diffusive motion, this expectation value can beused as a lower bound of the mean-squared displacement required to identify a particleas moving processively. In principle, a filament could exhibit (cid:104) ∆ x ( t ) (cid:105) (cid:54) = 0 due to theprocessive motion of other nearby filaments, but the results of such hydrodynamic effectsare unclear and they are ignored in the following derivations.Due to the anisotropy in the shape of the ellipsoidal myosin II filaments, it isexpected that they will exhibit different diffusivities in the directions parallel ( D a ) andperpendicular ( D b ) to their orientation (or semi-major axis). This means that, when inthe diffusive state, a filament will move with an average diffusivity D avg = ( D a + D b ) / x and y axes in the lab frame ( D x and D y respectively)will only tend to the value D avg at long times [9]. The analytical form of the diffusiontensor for an ellipse shows that D x and D y are functions of the initial orientation ofthe ellipse, and decay monotonically towards the average diffusivity D avg with the fficient tracking of myosin II filament position and orientation over time τ D = 14 D θ , (1)where the angular diffusivity D θ implicitly includes information about the averageeccentricity of the myosin II filaments. This is by definition the timescale for the diffusiontensor to become isotropic.Neglecting any single frame detection, the average dwell time of the myosin IIfilaments at the lowest ATP concentration (corresponding to the longest average dwelltime) was (cid:104) t (cid:105) = (1 . ± . D θ = (0 . ± . s − . This corresponds to an angular decorrelation time of τ D = (3 . ± .
10) s. The measured angular diffusivity was higher for the data set at thehighest ATP concentration, as on average fewer processive regions of track were detectedthat have high correlation in their orientations. In this case τ D = (1 . ± .
07) s, whichis more comparable to the average dwell time for the data set (cid:104) t (cid:105) = (1 . ± . τ D > (cid:104) t (cid:105) for all sample data sets, individual filament trajectories will maintaina good degree of correlation with their initial direction along their full length, and thediffusivities D x and D y along the lab frame axes for individual filament trajectorieswill be different from the long-term average diffusivity D avg . However, their average,( D x + D y ) /
2, provides an unbiased estimate of D avg even for individual trajectories [9].In this work, we sample a large population of trajectories whose initial orientations areisotropically distributed in the lab frame. In this case, even the single-axis diffusivities D x and D y are equal to D avg when the average over the whole set of trajectories isconsidered. The average diffusivity that we report, (cid:104) ( D x + D y ) / (cid:105) , is averaged over theensemble of all recorded trajectories, and is an unbiased estimate of D avg .Tracks with processive regions were isolated by requiring the mean-squareddisplacement of a myosin II filament to be greater than 16 D avg ∆ t (= 4 (cid:104) ∆ x ( t ) (cid:105) D avg )over a period of ∆ t = 2 s at some point along its track. This sets the requirement thata filament must have a dwell time of at least 2 s to be defined as moving processively,but as processive motion has been connected to stronger binding to the actin networkand longer dwell times [2] we believe this threshold to be reasonable. The value ofthe diffusivity D avg used to recognise processive tracks was derived from the sampledata set with the highest ATP concentration (using the same method as for the valueof D θ above and accounting for the 2D system), in order to minimise the effects ofprocessive motion on the measurement. It has been assumed that the average filamentdiffusivity is independent of ATP concentration. The effect of processivity on themeasured (translational) diffusivity and angular diffusivity for each of the sample datasets is shown in figure 3. The value of the diffusivity used to recognise processive trackswas D avg = (0 .
000 48 ± .
000 04) µ m s − , resulting in a threshold for processivity of16 D avg ∆ t ∼ . µ m . fficient tracking of myosin II filament position and orientation over time (cid:104) dr t +∆ t · dr t (cid:105) = 0 (where dr t = r ( t + ∆ t ) − r ( t )), and the variance in this quantity is,Var( dr t +∆ t · dr t ) = 8 D avg ∆ t , (2)after averaging over all possible initial orientations. This means that the standarddeviation in the scalar product dr t +∆ t · dr t is σ ( dr t +∆ t · dr t ) = 2 √ D avg ∆ t , even fora spatially anisotropic particle. We therefore require that the scalar product betweendisplacements at adjacent timesteps must be greater than the value of σ ( dr t +∆ t · dr t )for a purely diffusive particle in order for a point to be defined as processive.In order to minimise the number of false-positive detections of processive points ontracks, the smoothed, average signal, C + t = dr t +2∆ t · dr t +∆ t + dr t +∆ t · dr t , (3)has been calculated as in the work by Jeanneret et al. [8], and a myosin II filament isdefined as moving processively at time t if C + t > √ D avg ∆ t . Similarly, because thedecision of whether or not a filament is moving processively at time t should be invariant fficient tracking of myosin II filament position and orientation over time Figure 4: (a) Snapshots from iSCAT video with the lowest ATP concentration showingtrack of detected myosin II filament (cyan) with elliptical filament shape (magenta),labeled by time since the start of the video in seconds. (b) Track of example myosin IIfilament split into processive (blue) and diffusive (red) regions. Arrows show filamentorientation at each time. (c) Parameters used to define whether a point on the trackis processive or diffusive. Regions where either C + t or C − t (Eqs. (3 & 4)) are greaterthan 2 √ D avg ∆ t for at least 5 frames are processive, and b t = 1. (d) Distance frominitial position increases approximately linearly in processive regions, but plateausin long diffusive region. (e) Difference between orientation angle and direction ofpropagation (calculated from Eq. (7)) appears to have more frequent fluctuations withlarger amplitudes in the diffusive regions.under time reversal symmetry, for each point on the track the averaged correlation C t iscalculated for both directions in time. Time inversion generates a new form of Eq. (3) fficient tracking of myosin II filament position and orientation over time C − t = dr t − t · dr t − ∆ t + dr t − ∆ t · dr t , (4)but the value to compare this to for a purely diffusive particle is invariant under thistransformation, σ ( dr t − ∆ t · dr t ) = 2 √ D avg ∆ t . An example of a track that has beensplit into diffusive and processive regions is shown in figure 4.Using the average correlation in both directions in time to define whether a myosinII filament is moving processively ensures that points on a filament’s track at the edgeof a processive region are not neglected incorrectly (for example C + t cannot be definedfor the final three points of any track), and that short ( ≤ b t hasbeen defined such that b t = 1 at time t if either C ± t > √ D avg ∆ t , and b t = 0 otherwise,as shown in figure 4(c). The SEP Python package generates the orientation angles of each elliptical particle itdetects in an image (above the threshold area and intensity described in section 2.2)on a − π/ ≤ θ ≤ π/ − π ≤ θ ≤ π domain instead. Extending this domain further and tracking filamentorientation on the domain −∞ < θ < ∞ ensures that no large jumps are observed infilament orientation due to the periodic boundaries of a − π ≤ θ ≤ π domain. Thisre-parameterisation is especially important for the calculation of angular diffusivity.Of the two possible initial orientation angles (parallel and anti-parallel to thesemi-major axis of the elliptical myosin II filament detection), it is assumed thatthe correct choice will minimise the average (over time) of the difference betweenthe filament’s orientation at time t , θ t , and propagation direction at time t , φ t =arctan((¯ y t +∆ t − ¯ y t ) / (¯ x t +∆ t − ¯ x t )) (where ¯ x t , ¯ y t are the co-ordinates of the centroid ofthe detected ellipse at time t ). This average is calculated over either the processiveregion(s) of track, or the entire track if the filament exhibits purely diffusive motion. Ithas been found that filaments moving processively have a smaller average value of φ t − ¯ θ t fficient tracking of myosin II filament position and orientation over time φ is aligned with orientation angle θ + for majority of track.Black data shows the orientation extracted by the SEP Python package, red and bluedata represent possible orientation tracks in angle-space, θ + and θ − respectively. The θ + data maximises the sum in Eq. (5). (b) Purely diffusive track, propagation direction φ varies faster and with greater amplitude than orientation angle θ ± . Orientationsextracted by the SEP Python package are split amongst both of the possible orientationtracks in angle-space.(where ¯ θ t = ( θ t +∆ t + θ t ) /
2) [2], with examples shown in figure 5, so using processiveregion(s) to derive the filament’s initial orientation is preferred if possible.Once a detection is added to a track it is assigned the orientation angle (orientedeither parallel or anti-parallel to its semi-major axis) that minimises the filament’schange in orientation between frames. It has been assumed that a filament cannot rotateby more than ∆ θ max = π/ −∞ < θ < ∞ domain. It has also been assumed thatthe propagation direction of a filament cannot vary by more than ∆ φ max = π betweenframes, so the new propagation direction minimises the change in the directions betweenframes taking into account the 2 π periodicity of the domain.Of the two possible orientation tracks in angle-space (that result from the twoavailable initial filament orientations) the selected path is the one that maximises thesum, Ψ( θ, φ, T ) = ( T/ ∆ t ) − (cid:88) i =0 cos( φ i ∆ t − ¯ θ i ∆ t ) , (5)where T is the total duration of the iSCAT video being analysed (with correspondingtime between frames ∆ t ). If the filament exhibits processive motion, then the sum inEq. (5) is instead taken from i = T s / ∆ t to i = ( T f / ∆ t ) −
1, where T s and T f are the startand finish times of the processive region(s). The cos term inside the sum is a maximumwhen the average filament orientation is aligned with its propagation direction, and is a fficient tracking of myosin II filament position and orientation over time π periodicity of the domain. The values Ψ ± ( θ ± t , φ t , T ) corresponding to the two possibleorientation paths θ ± t (labelled as one path will always start with a positive orientationangle and the other a negative one) will always be separated by a factor of −
1, as,Ψ ± ( θ ± , φ, T ) = cos( ± π )Ψ ∓ ( θ ∓ , φ, T ) = − Ψ ∓ ( θ ∓ , φ, T ) , (6)so only one path needs to be followed in order to choose the correct initial orientation.By exploiting the domain of the arccos function from the Math Python package,the magnitude of the difference between the orientation angle and the direction ofpropagation can be calculated at each time t as,Φ(¯ θ t , φ t ) = | arccos(cos( φ t − ¯ θ t )) | . (7)Using Eq. (7), the distance moved by a myosin II filament at each time, | ∆ x t | = r t , can be separated into components parallel, | ∆ x t | para = r t | cos(Φ(¯ θ t , φ t )) | , andperpendicular, | ∆ x t | perp = r t | sin(Φ(¯ θ t , φ t )) | , to the filament orientation. The error inindividual orientation measurements due to image pixelisation and inherent pixel noisecan be estimated by using the tracking procedure presented here to observe particlespermanently stuck to a surface.Using the results of this angle tracking, a histogram of ¯ θ t − ¯ θ can be plotted foreach time t after the initial observation of a filament. These results show that theorientation distribution evolves as a Gaussian with width proportional to time ∼ D θ t (results not shown), and hence that it was correct to use a linear fit to calculate theangular diffusivity from the mean-squared angular displacement data.
3. Discussion
The computational methods described here were developed with the aim of analysingthe motion of myosin II filaments containing multiple motor protein domains whenbound to an underlying actin filament network. Regular point or circular body trackingalgorithms could not be used to accurately parameterise the motion of the elongated,ellipsoidal myosin II filaments, and would not be able to provide any information abouttheir angular fluctuations. Splitting tracks into regions of diffusive or processive motionallows us to potentially probe the dynamics of different bound states for the myosin IIfilaments, which could be extended to the study of different biological systems. Followingthe development and calibration of the SPT method outlined in this paper, analysis ofthe full myosin II filament data set has yielded interesting information about the dwelltime distribution and spatial dynamics of bound myosin II filaments as a function ofthe ATP concentration of the system [2]. This novel SPT method can be used to studythe dynamics of particles at varying length-scales and has potential applications in thefields of fluorescence and light microscopy. Tracking and analysis code is available uponrequest. fficient tracking of myosin II filament position and orientation over time Acknowledgments
DVK acknowledges the Kukura lab for enabling the imaging of myosin II filamentdynamics on their iSCAT microscopes. DVK thanks the Warwick-Wellcome QBP forfunding. LM and MP gratefully acknowledge support from Leverhulme Trust GrantRPG-2016-260.