Automated visual tracking for studying the ontogeny of zebrafish swimming

Capturing the detailed locomotory kinematics and behavior of biological organisms plays an important role in a wide variety of research disciplines(Fry et al., 2003; Gray, 1933; Park et al., 2001; Pfau et al., 2005). Many studies in biomechanics, neuroethology and developmental biology rely on the analysis of video sequences. Such analyses stand to benefit greatly from computer vision techniques, which can accurately and efficiently screen complex morphometric and locomotory characteristics.

Following the pioneering work of Eadweard Muybridge(Muybridge, 1887), scientists have developed a wide range of methods to record and quantify animal locomotion. Many studies have relied on tracking landmarks, such as leg joints(Full and Tu, 1991) or markers on the animal's body (Hedrick et al.,2004; Standen and Lauder,2005). However, it is not always feasible or advisable to attach markers, e.g. during field observations or on animals that have a mucous layer(e.g. many fish) or motile skin (e.g. horses). In these cases, tracking the silhouette of the animal or a body part can offer an alternative [e.g. silhouettes of ascidian larvae and swimming fish(McHenry, 2001; Tytell and Lauder, 2002)].

Tracking swimming movements by filming the swimmer's silhouette has a long tradition (Batty, 1984; Gray, 1933). Silhouettes lend themselves to automatic tracking (Tytell and Lauder, 2002). Furthermore, tracking routines can be refined by using information about the shape of the animal or limb to be tracked [e.g. insect wings (Fry et al.,2003; Willmott and Ellington,1997)]. Error reduction as a result of including shape information is particularly valuable when studying the propulsive body wave of fishes. The exact shape of a fish's body wave depends on the interaction between the solid body of the fish and the surrounding water. Exact knowledge of the wave shape can provide valuable insight into swimming performance and the interaction between undulating body and surrounding water(Cheng and Pedley, 1998; McHenry et al., 1995; Wainwright, 1983). For example, the stiffness of the fish's body affects the shape of the body wave and thereby influences the propulsive performance of the fish. An interesting case study is presented by the stocksteif mutation in zebrafish,which is characterized by an overossification of the notochord. The axial skeleton of these mutants is a stiff, bony rod, contrasting the flexible series of articulating vertebrae in their wild-type siblings.

The zebrafish has long served as a convenient model to study the various aspects of fish swimming (Fuiman and Webb,1988; Thorsen et al.,2004). Automated tracking and analysis systems have been previously developed for zebrafish (Bang et al., 2002; Blaser and Gerlai,2006). However, these systems track the fish only as a point and cannot quantify body wave kinematics of swimming. Other studies of zebrafish swimming have manually tracked the fish to quantify the body posture, a method that is both time-consuming and potentially prone to subjective errors(Budick and O'Malley, 2000; McElligott and O'Malley, 2005; Müller and van Leeuwen,2004). Tytell and Lauder used a semi-automated method to estimate the fish midline by manually indentifying the snout and tail and automatically estimating the midline from the extracted silhouette(Tytell and Lauder, 2002). Other authors have relied on `skeletonizing' algorithms that dissolve a binary image representing the animal's silhouette down to its midline(Cronin et al., 2005; Geng et al., 2004; McHenry, 2001). Although this approach is automated, it will not estimate the correct midline if other objects are present in the binary image because it cannot distinguish between pixels that belong to the animal's silhouette and those that belong to a different object. As a result, these algorithms will not correctly estimate the fish's body posture when there is environmental clutter such as other fish, plants, or a hair used to initiate behavioral responses, as we did in our recordings. Automated kinematic analysis of multiple zebrafish larvae was recently demonstrated. However, this particluar analysis technique utilizes an image filter that is customized for the appearance of zebrafish larvae of a specific age (Burgess and Granato,2007). This technique does not extend nicely to zebrafish of different ages, other fish species, or when environmental clutter is present.

Here, we present a complete method for accurately and efficiently quantifying the body posture of zebrafish and other organisms with symmetric medial profiles. Our approach directly models the shape of the animal and utilizes locations of high contrast in the image to estimate its posture. The posture estimate is calculated using techniques that remain robust to clutter. The detailed swimming motion is estimated based on dorsal images of the fish recorded at sufficiently high frame rate(Harper and Blake, 1989),which enables a quantitative evaluation of the animal's kinematics and dynamics (provided the mass distribution of the animal is known). In the first section, we review the overall approach for performing model-based visual tracking. Subsequently, we develop a detailed geometric model for the zebrafish, including appropriate motion and measurement models that use information from the previous and current frame to estimate the fish's current position and posture. Finally, we demonstrate the capabilities of our tracking approach on zebrafish performing an escape response at three stages during their development (from larvae to juvenile).

Zebrafish (Danio rerio Hamilton 1822) eggs were collected after mating one stocksteif heterozygous female with one stocksteif heterozygous male. The batch of eggs contained both stocksteif mutant and wild-type embryos, but the mutant phenotype does not become apparent until 5 days post fertilization (d.p.f.). The embryos were reared at the optimal rearing temperature of 28°C. After hatching,the embryos were fed Paramecium (5 and 6 d.p.f.) and Artemia(from day 7 onwards). A fast startle response was recorded at 5, 15 and 28 d.p.f. for both wild-type and stocksteif mutant animals using a high-speed video camera (Photron, APX RS, 1500 frames s^–1,1024×1024 pixels, exposure time 1/8000 s) fitted with a 105 mm Nikon lens. The startle responses were elicited by touching the animals with a horse hair. We analyzed recorded sequences from the initiation of the escape response to the moment when the fish either leaves the field of view or ceases active swimming. Our sequences therefore include stages 1 and 2 of an escape response, and usually several tail beats that are part of stage 3(Weihs, 1973) (Movie 1 in supplementary material).

Recently, the Kalman filter has been considerably improved by statistical linearization and now performs well when applied to nonlinear motion and measurement models (Ito and Xiong,2000; van der Merwe and Wan,2003; Nörgaard et al.,2000; Sibley et al.,2006). These Sigma Point Kalman Filters (SPKF) recursively estimate the optimal fish parameters and their uncertainty by taking the weighted average of a set of regression points drawn from the probability density function that describes the fish's shape and location at the previous time step. In this study, we leverage the power of SPKFs to solve our nonlinear estimation problem. However, in order to estimate the fish parameters, we must also develop equations for the motion model f and the measurement model h that make sense for a swimming zebrafish during development. In the following sections we apply a two-step approach. First we construct a flexible geometric model of a zebrafish outline. Second,we develop appropriate motion and measurement equations for this model, which are used in conjunction with the SPKF to automatically optimize the fish parameters based on the processed image data.

Our complete fish model is a symmetric planar mesh with grid parameters u and v, which define a local coordinate frame in the fish body in the length and width directions, respectively. The points in this mesh[x(u,v) y(u,v)]^T are continuous functions of u and v, so we can evaluate them at whatever discrete grid we choose. For our experiments, we sample 30 uniformly spaced points in the domain [–(1–γ)L,γ L] for u and three uniformly spaced points in the domain [–1, 1] for v. Because the number of sampled points in the (u,v) grid remains constant during tracking, the mesh points that make up our fish model become functions of the position and shape parameters, p. We construct the centerline of our symmetric mesh by integrating the unit tangent vector

Creating deformable models based on medial profiles has been used in segmentation problems in medical imaging(Hamarneh and McInerney, 2001)and for tracking multiple C. elegans from microscopy images(Fontaine et al., 2006; Roussel et al., 2007). Because zebrafish are laterally symmetric about their body axis, our medial profile representation offers several advantages for the tracking framework. Each B-spline basis function is defined only over a subregion of the fish body. Therefore, the local bend amplitudes, α_j, have local control over the degree of bending in the body. This property is analogous to the fish anatomy, where contractions of individual muscles affect the bending over subregions of the fish's body. In summary, our parameterization of the centerline has few degrees of freedom, requires no training data, and offers a natural and anatomically sound way to constrain the fish's length and designate certain regions as stiff.

The motion model of the fish predicts the fish parameters at the current time step based on the parameters calculated at the previous time step. By predicting the motion, it provides a better initial estimate of the fish parameters before they are updated using the measurement model. We assume that the fish movement is a combination of undulatory motion along its body and a displacement of the whole body along the centerline axis. The undulatory motion of steady swimming in fish consists of a traveling wave of increasing amplitude from head to tail. Given our model of the zebrafish, this motion is governed by the change in local bend amplitudes, α_j,from the previous time step to the current one. Given that our video sequences are filmed at a sufficiently high frame rate, we make the assumption that the current local bend amplitudes are equal to the previous ones plus a random variable drawn from a zero mean normal distribution. This is a simple way to allow in the model that we are more uncertain about the body wave shape of the fish since it has moved between frames. We further assume that the fish has constant velocity, η₁, between frames during displacement along its body axis that is corrupted by acceleration noise, η₂. This velocity is also estimated and thus incorporated into the set of fish parameters

The present tracker system assumes that we are able to measure the location of the fish boundary, so our measurement model consists of the boundary points, q_i=[x_iy_i]^T,along with their outward normal vectors,

Any tracking algorithm relies on an initial estimate of the object location. To achieve this, we have developed a semi-automated initialization routine that operates on the first movie frame and extracts three important pieces of information for tracking: (1) an estimate of the background image,(2) the initial fish parameters p₀, and (3) the width function R(u). The background model is estimated by selecting a region around the fish and erasing it using the built-in Matlab™ function `roifill', which smoothly interpolates inward from the pixel values on the boundary of the user-defined region. Next, the background image is used to segment the first movie frame, and the Matlab™ function `bwboundaries' calculates the fish boundary from the resulting binary image. The user is then requested to click on the snout and tail locations of the fish, which allows us to divide the boundary into the left and right discrete boundaries of the fish denoted B_L(i) and B_R(i),respectively.

The initial centerline of the fish will be the curve that is equidistant from the left and right boundaries. To determine this centerline, we use a modification of the integral area distance used in(Roussel et al., 2007) which iterates and converges to the fish centerline. Our method works by first finding the closest point in B_L to each point in B_R and vice versa. Then, we calculate the median locations between each set of corresponding points in B_Land B_R. These median locations are assigned as the new left and right boundaries and this process is repeated until the boundaries converge onto the true discrete centerline C(i) (see Appendix for more details). The estimated centerline curve C(u) is determined by fitting a B-spline through the points calculated from C(i), and is illustrated in Fig. 6A along with the original left and right boundaries it was derived from. The bend angle function

In this section we present automated tracking results for wild-type and stocksteif mutant zebrafish at 5, 15 and 28 d.p.f. Fig. 7 shows the raw centerlines of the zebrafish estimated by the tracker at fixed time intervals and demonstrates the quality of the proposed method. Our method successfully tracks fast escape responses of fish larvae(Fig. 7A,D) despite occasional partial occlusions [in this case, the hair used to induce the escape response(Fig. 7E)]. Furthermore,tracking takes an average computation time of 5.5±1.7 s frame^–1 on a 3.0 GHz Intel® Xeon processor, which enables us to analyze large datasets very quickly. To calculate an upper bound on the accuracy of our tracking approach, we created a synthetic movie sequence by rendering our model along a known trajectory. We then tracked this synthetic sequence and measured the average error between the known and estimated centerlines over time (Fig. 8). We were able to localize the synthetic data within 0.5% of the body length on average; our errors in estimating the real movie sequences will be slightly larger because they contain additional noise sources not present in the synthetic one.

From these centerlines, we wish to extract important kinematic parameters to gain insight into developmental influences on the propulsion mechanisms of swimming fish, and vice versa: the mechanical influences on the development of the fish. One of these parameters is body axis curvature, which provides information about the muscle strains the fish undergoes. To measure curvature, we apply spatial smoothing to the extracted centerlines and then apply temporal smoothing directly to the curvature values. Spatial smoothing is performed by fitting a cubic B-spline curve to the extracted centerlines,and the curvature is calculated directly from these smoothed curves (see Appendix for details). We calculate curvature from these smoothed centerlines instead of deriving it from the raw centerline of our geometric model because the model contains a discontinuity in the curvature at the location where it becomes stiff.

Temporal smoothing is performed by applying a low-pass filter to the curvature values across time at fixed locations along the fish's body (we use 51 uniformly spaced points). The cut-off frequency for the filter is chosen based on visual inspection of the magnitude response of the curvature's Fourier transform at each of the body locations. Fig. 9 plots the error between the filtered and unfiltered centerlines and curvature values for both wild-type and stocksteif zebrafish. With average errors around 0.1%of body length and 0.1 for the centerline and curvature, respectively, our filtering process retains all-important information and does not compromise the accuracy achieved by the tracker. The resulting curvature profiles after filtering are shown in Fig. 10. Comparing stocksteif with wild type at age 15 and 28 days (Fig. 10B,E, Fig. 10C,F) we find that the peak curvature values of stocksteif are smaller than those of wild type.

To analyze the performance of our fish, we measure several additional kinematic parameters. The angular acceleration is the second temporal derivative of the fish bend angle function,

Typical measures of escape swimming performance include displacement,speed, and acceleration of the fish center of mass (COM) when stretched straight (Domenici and Blake,1997; Walker,1998). To estimate the location of the COM, we reconstructed the fish volume from dorsal and lateral photographs assuming an elliptic cross-section (McHenry, 2001)(see supplementary material Figs S1–S3). The center of volume (COV) and center of area from the dorsal view (COA) were calculated. The COV will correspond to the COM assuming the fish has uniform density. COV calculations were only performed on wildtype zebrafish for comparison with published literature on escape responses. However, COV calculations do not lend themselves to high-throughput analysis because separate photographs of each individual fish must be acquired to account for the variation in morphometry for fish at a given age (very large variation exists for the stocksteif at a specific age). Instead, we propose to use the COA as a location for comparison between wild type and stocksteif because it is easily measured from the video sequence and has similar speed and acceleration profiles to the COV (see Fig. 13). Fig. S4 in supplementary material demonstrates that the COA location has a maximum deviation of 5–6% body length from the COV at age 5 days and the deviation decreases as the fish gets older. Fig. 13 illustrates the results of these measurements. Speed and acceleration are calculated using the MSE spline method (Walker,1998). The stocksteif consistently exhibits lower peak acceleration and speed compared with wild type at each age. Although these preliminary measurements indicate consistent discrepancies between the wild type and stocksteif due to the stiffer vertebral column present in the mutant, a larger sample must be analyzed to produce statistically significant results, such an analysis will be presented in our next paper.

We have presented a method for estimating the body posture of zebrafish within laboratory environments using flexible geometric models and nonlinear estimation. Given the generalized mathematical framework used to model the fish's appearance, this method should track any fish species with a symmetric medial profile and that swims by undulating the body. The discrete time,dynamic state, space model also provides a general framework for performing statistical inference that is robust to outliers and enables tracking during partial occlusions. This is a strong improvement over previously developed methods that are either manual (Budick and O'Malley, 2000; McElligott and O'Malley, 2005; Müller and van Leeuwen, 2004), require a perfectly segmented image with no environmental clutter (Cronin et al.,2005; Geng et al.,2004; McHenry,2001), or are customized for fish of a specific size and appearance (Burgess and Granato,2007).

The assumption of planar motion is a limitation in the proposed method,which arises directly from the recorded material – top-view shots with a single camera. For behaviors that contain large out-of-plane motions, the geometric model will not accurately represent the appearance of the fish. For such behaviors, a 3D version of our model-based tracker would be required; our principle is sound, but should be extended. The tracker may also fail if a significant portion of the fish becomes completely occluded by environmental clutter (e.g. the entire head). Nevertheless, this could be improved by extending the algorithm to use a more advanced nonlinear estimator, a more advanced motion model, or a direct model of the occlusions. However, in laboratory settings, where many environmental parameters are controlled, this technique represents an accurate and fully automatic approach to quantify behavior and will facilitate studies requiring the analysis of many and long image sequences.

In addition to traditional swimming performance indicators we explored two new ways of analyzing the kinematics data: by plotting the angular acceleration as a function of time and the frequency response along the body. We also used more objective mathematical definitions and corresponding algorithms to quantify standard variables such as tail beat frequency, wave speed and length. The center of area (COA) of the fish dorsal view was proposed as a valid location for comparison between wild-type and stocksteif fish because its ease of measurement lends itself to high-throughput analysis. Our preliminary comparison between wild-type and stocksteif swimming performance indicators already suggest significant differences. Hence, this preliminary data analysis of swimming fish illustrates the capabilities of the automatic fish tracker and bodes well for gaining a complete understanding of how stiffness of the vertebral column affects swimming performance.

LIST OF ABBREVIATIONS AND SYMBOLS

 
• B_L,R

  initial left and right boundaries of fish image (in pixels)

 
• C(i)

  initial centerline of fish image (in pixels)

 
• C(u)

  initial centerline of fish image modeled as B-spline

 
• COA

  center of area

 
• COM

  center of mass

 
• COV

  center of volume

 
• \({\vec{e}}_{1}\)

  tangent vector to fish model centerline

 
• \({\vec{e}}_{2}\)

  normal vector to fish model centreline

 
• f(.,.)

  motion model of zebrafish

 
• h(.,.)

  measurement model of zebrafish

 
• H(p)

  complete geometric model of zebrafish

 
• L

  total length

 
• n_i

  outward normal vectors of model boundary

 
• \(N_{{\Theta}}\)

  number of B-spline basis functions used to model centerline bend angle

 
• p₀

  initial fish parameters

 
• p_k

  fish parameters at the kth time step

 
• P_ab

  covariance matrix associated with the random variables a and b

 
• q_i

  model boundary point

 
• r_i

  edge feature point in image

 
• R(u)

  width profile of zebrafish model

 
• \({\vec{T}}\)

  global translation vector of fish model

 
• (u,v)

  grid parameterization of model mesh

 
• z_k

  image measurements at the kth time step

 
• αj

  local bend amplitudes

 
• \({\vec{{\alpha}}}\)

  bend angle amplitude

 
• β

  model scaling term (pixels mm^–1)

 
• γ

  ratio of stiff region of zebrafish model to entire length

 
• Δt

  time between frames

 
• η0

  axial displacement of fish between frames

 
• η1

  axial velocity of fish between frames

 
• η2

  axial acceleration fish between frames, sampled from normal distribution

 
• \({\Theta}\)

  (u)

  bend angle of centerline in zebrafish model

 
• λ_A(B)

  nearest-neighbor function; the element in A that is closest to B

 
• Λ(u)

  s by

  \(N_{{\Theta}}\)

  matrix of B-spline basis functions evaluated at s sampled grid points in u

 
• \({\Phi}_{j}^{k}\)

  (u)

  B-spline basis function of order k

The authors would like to thank B. Walderich, H. M. Maischein and J. Odenthal of the Hubrecht laboratory for graciously allowing us to use video of the stocksteif mutant for demonstrating the ability of our tracker system. We also thank Ansa Wasim for recording the movie sequences and Henk Schipper for capturing the photographs used in the COV calculations. E.I.F.,A.H.B. and J.W.B. thank the Beckman Institute of Caltech for providing financial support of this project. D.L. and J.L.v.L. are funded by NOW-ALW grant 817.02.012.