Abstract
Stereotyped behaviors are series of postures that show very little variability between repeats. They have been used to classify the dynamics of individuals, groups and species without reference to the lowerlevel mechanisms that drive them. Stereotypes are easily identified in animals due to strong constraints on the number, shape, and relative positions of anatomical features, such as limbs, that may be used as landmarks for posture identification. In contrast, the identification of stereotypes in single cells poses a significant challenge as the cell lacks these landmark features, and finding constraints on cell shape is a nontrivial task. Here, we use the maximum caliber variational method to build a minimal model of cell behavior during migration. Without reference to biochemical details, we are able to make behavioral predictions over timescales of minutes using only changes in cell shape over timescales of seconds. We use drug treatment and genetics to demonstrate that maximum caliber descriptors can discriminate between healthy and aberrant migration, thereby showing potential applications for maximum caliber methods in automated disease screening, for example in the identification of behaviors associated with cancer metastasis.
Introduction
Moving in the right way at the right time can be a matter of life and death. Whether avoiding a predator or searching for food, choosing the correct movements in response to specific stimuli is a crucial part of how an organism interacts with its environment. The repetitive, highly coordinated movements that make up behavioral stereotypes have been shown to be entwined with survival strategies in a number of species, for example the incredible correlation in posture between prey capture events in raptors^{1} and the escape response of C. elegans when exposed to extreme heat^{2}. Understanding these stereotypes is vital to creating a full picture of a species’ interactions with its environment. If stereotypes represent evolved, selectiondriven behavior in animals, might the same not be true for singlecelled organisms?
This point of view may be particularly useful in understanding chemotaxis, the guided movement of a cell in response to a chemical gradient. During chemotaxis, eukaryotic cells change their shape through the repeated splitting and extension of actinrich structures called pseudods^{3,4,5}. Though this behavior is well known, the study of chemotaxis has traditionally focused on the signaling events that regulate cytoskeletal remodeling. Even where pseudopods are acknowledged to be relevant, the focus is on the biochemical mechanisms that generate and regulate them^{6,7,8}. These mechanisms are, however, staggeringly complex^{9} and the way chemotaxis emerges from these lowerlevel processes remains largely unknown. Rather than delving deeper into the network of biochemical interactions, we can instead learn from the shape changes and movements that this intricate machine has evolved to produce. Such an approach, also known as morphological profiling, shows great promise in biomedicine^{10}.
Here, we explore this question using Dictyostelium discoideum, a model for chemotaxis noted for its reproducible directional migration towards cyclic adenosine monophosphate (cAMP)^{11,12}, which it senses using typical Gprotein coupled receptors. To capture cell shape (or posture) at any given point in time, we employ Fourier shape descriptors, a rotationally invariant method of quantifying shapes used previously to show that cell shape and environment are intrinsically linked^{13} (Fig. 1A). These shape data are naturally of high dimensionality, making further analysis difficult. We reduce their dimensionality using principal component analysis (PCA), a method used previously to obtain the key directions of variability from the shapes of cells^{13,14,15} and animals (Fig. 1B)^{2,16,17}. Our final challenge (and the focus of this paper) is to quantify behavior, which we define as the movement between shapes. There are many potential ways to do so^{18,19}, however we have adapted the variational maximum caliber (MaxCal) approach^{20,21} to this end. These methods have several advantages over conventional alternatives: Firstly, Fourier descriptors capture all available information on shape, and the subsequent PCA provides a natural quantitative means of discarding its less informative elements. Easier methods, such as measuring aspect ratio or eccentricity, require us to assume that they provide a useful description a priori, and cannot inform us how much (or what) we have discarded for the sake of simplicity. Secondly, our chosen methods are blind to the researcher’s preconceptions, as well as to previous descriptions of shape and behavior. Any behavioral modes identified have emerged from the data without our deciding on (and imposing) them, as we might if using supervised machine learning or fitting parameters to a preconceived biochemical model. Finally, the minimal model we construct using maximum caliber makes no reference to any underlying biochemistry and therefore cannot make potentially incorrect assumptions about it. We demonstrate the usefulness of these methods by showing that they successfully discriminate between the behavior of drugtreated or genetically altered cells and their parental strains.
Results
Maximum caliber approach to behavioral classification
Cells continuously change shape as they migrate, creating trajectories in the space of shapes that are specific to their situation. For example, we have previously shown that cells follow different shape trajectories in environments with low and high chemoattractant signaltonoise ratios^{13}, here defined as the local gradient squared over the background concentration (Fig. 1C). In this example, it is important to note that the distributions of cell shape for each condition overlap significantly. This means that it is not always possible to accurately determine the cell’s condition from a static snapshot of its shape. In contrast, the dynamics of shape change in each condition are clearly distinct. Our aim here is to quantify the details of these shape changes, making a small set of values that can act as a signature for a given mode of behavior. We can then use such signatures to quantitatively compare, or to discriminate between, various conditions or genotypes. To this end, we employ the MaxCal method (Fig. 2A).
MaxCal was originally proposed by E. T. Jaynes^{22} to treat dynamical problems in physics, and is much like his betterknown maximum entropy (MaxEnt) method used in equilibrium physics. The motivation is the same for both; we wish to find the probability of a given structure for a system in a way that neither assumes something we do not know, nor contradicts something we do know, i.e. an observation of the system’s behavior. In the case of MaxEnt, this is achieved by finding the maximum Shannon entropy with respect to the probabilities of the set of possible configurations the system can take^{23}. MaxCal uses the probabilities of the possible trajectories a dynamical system can follow instead. In this case, the entropy is replaced by the caliber C^{20}, so called because the flow rate in a pipe relates to its caliber, or internal diameter. In essence, the method extracts the degree to which different rates or events within the system are coupled, or cooccur beyond random chance. This method has previously been used to sucessfully predict the dynamics of neuronal spiking, the flocking behavior of birds and gene regulatory networks^{24,25,26}.
Generally, the caliber takes the form
where p_{j} is the (potentially timedependent) probability of the jth trajectory. The first term on the righthand side of Eq. (1) represents a Shannonentropy like quantity and the second ensures that the p_{j} are normalized. The third constrains the average values of some properties S_{n,j} of the trajectories j to the values of some macroscopically observed quantities 〈S_{n}〉, making sure we do not contradict any known information.
By maximizing the caliber, the probabilities of the trajectories
are found, where \(Q={\sum }_{j}\exp (\sum _{n}{\lambda }_{n}{S}_{n,j})\) is the dynamical partition function and {λ_{n}} are a set of Lagrange multipliers. This Boltzmanntype distribution fulfils detailed balance, even for nonequilibrium problems. Practically, the problem is to find these Lagrange multipliers (and hence, the partition function). To this end, we exploit their relations to the externally observed average values of some quantities
where the values of the 〈S_{n}〉 are determined from experiment. This training process is equivalent to maximumlikelihood fitting to observed data.
As our interest is in cell shape and motility, we derive our values for the 〈S_{n}〉 from the shape dynamics of migrating Dictyostelium cells. In order to effectively parameterise our model, we must constrain the continuum of possible shape changes to a much smaller set of discrete unit changes in our principal components (PCs). We therefore build our model from discretized values of the shape measures PC 1 and 2, assigning them to the variables N_{1} and N_{2}, respectively. Their values are analagous to particle numbers in a chemical reaction. The switching between continuous and discrete variables is possible as \(\frac{{\sigma }_{x}}{\langle {N}_{x}\rangle }\approx 0.035\) is small with x = 1, 2 for PC_{x} and σ_{x} the standard deviation. We reduce the size of the timestep δt until we no longer observe changes greater than 1 in a single δt (similar to deviations of master equations). As PC 2 accounts for less overall variation than PC 1, (see Fig. 1B), we naturally reach this minimal change for a much larger value of δt, which is undesirable because by the time δt is small enough for PC 1, changes in PC 2 are almost never observed, making correlations between the two PCs difficult to detect. We therefore scale all changes in PC 2 by a factor of σ_{1}/σ_{2} in order that unit changes are observed in both PCs for a similar value of δt. Practically, our training data yielded a δt of 0.1875 s (as each 3 s frame in the video data was divided into 16 sections, in which the intermediate PC values were linearly interpolated).
The advantage of limiting the possible macroscopic shape changes in δt to the following: an increase, a decrease, or no change in each PC. As changes in each PC can be considered independently, this gives us a total of 3 × 3 = 9 cases (that is, no change form the current position, or a move to any of the 8 neighbouring spaces, see Fig. 2A inset). These macroscopic cases are taken to be the observable effects of an underlying microscopic structure. From our analogy of a chemical reaction, we treat increases to be particle creation and denote the microscopic variable for an increase in trajectory j as \({S}_{+,j}^{x}\), where x ∈ {1, 2} correspond to PC 1 and 2, respectively. For small δt this variable is binary, taking the value 1 when N_{x} increases over a single timestep and taking the value 0 otherwise. Decreases will be treated as particle decay, with N_{x} separate variables \(\{{S}_{,j}^{x,i}\}\) are used to denote decays for the ith particle, with 1 ≤ i ≤ N_{x}. These \(\{{S}_{,j}^{x,i}\}\) are equal to 1 if the ith particle decays in δt and equal to 0 otherwise. Hence, in each δt there are N_{x} + 2 possible microtrajectories for each component; an increase, no change, or the removal of any particle N_{x} (Fig. 2B). We choose such a firstorder decay over a zerothorder decay in order to introduce a virtual force, bringing the system back toward the mean (see Fig. S2). As the two components may change independently, there are (N_{1} + 2)(N_{2} + 2) possible microtrajectories in a single δt over PC 1 and 2. Applying Eq. 3, we constrain the probabilities of these microtrajectories such that they agree with the macroscopically observed rates \(\langle {S}_{\alpha }^{x}\rangle \), with α ∈ {+, −} an increase or decrease in component x, respectively.
We then expand the model to include a following timestep, allowing us to capture shortterm correlations between events. This increases the number of possible trajectories substantially. The number of microtrajectories in a given timestep depends on N_{x} at time t + δt, and this quantity is different dependent on the pathway taken in the first timestep, so we must include this history dependence. For example, a reduction in component x can happen in N_{x} ways, and will cause N_{x} to go down by one. This change is followed by (N_{x} − 1) + 2 possible microtrajectories in the following timestep. Multiplying the quantities for each timestep gives us N_{x}(N_{x} + 1) microtrajectories in which there is a decrease in the first timestep. Accounting for the effect of the changing values of N_{1} and N_{2} over interval t, t + δt in each microtrajectory on the interval starting at time t + δt, the number of microtrajectories over 2δt is \(({N}_{1}^{2}+3{N}_{1}+5)({N}_{2}^{2}+3{N}_{2}+5)\). Each observable has a corresponding value in trajectory j of \({S}_{\alpha \beta ,j}^{xy}\), which is 1 if the correlation is observed and 0 otherwise. We can reduce this to 10 timecorrelated observables by assuming symmetry under orderreversal, i.e. \({S}_{\alpha \beta ,j}^{xy}\equiv {S}_{\beta \alpha ,j}^{yx}\) (Fig. 2C). This assumption is justified: if we consider a negatively correlated movement between PC1 and PC2, we may see transitions in the order 1+, 2−, 1+. Here the two couplets 1+, 2− and 2−, 1+ both represent the same phenomenon (see Fig. S3). This leads to an additional 16 observables \(\langle {S}_{xy}^{\alpha \beta }\rangle \), where x, y ∈ {1, 2} are shape PCs and α, β ∈ {−, +} denote a change in the component displayed above. We constrain our analysis to the first two shape components only, as further components account for relatively little residual variance in shape, whilst increasing computational complexity geometrically.
As an example, we show the partition function in a single shape component, in which there are 5 observables, \(\{\langle {S}^{+}\rangle ,\langle {S}^{++}\rangle ,\langle {S}^{+}\rangle ,\langle {S}^{}\rangle ,\langle {S}^{}\rangle \}\):
where \({\gamma }^{\alpha }={e}^{{\lambda }^{\alpha }}\) corresponds to a rate (when divided by δt), with λ^{α} the Lagrange multiplier associated with observable 〈S^{α}〉. The first line in Eq. (4) shows all possible transitions that begin with an increase over the first timestep, and so the whole line shares the factor γ^{+}, the rate of increase. A subsequent increase contributes a further γ^{+}, as well as a coupling term γ^{++} which allows us to capture the likelihood of adjacent transitions beyond the naive probability γ^{+}γ^{+}. A subsequent decrease can happen in N + 1 ways, each linked to the rate of decrease γ^{−}. The term γ^{+−} is a coupling constant, controling the likelihood of an adjacent increase and decrease beyond the naive probability γ^{+}γ^{−}. Finally, the +1 allows for the possibility of no transition occurring in the subsequent timestep. The second and third lines correspond to a decrease in the first timestep, and no transition occurring in the first timestep, respectively.
The Lagrange multipliers corresponding to observables are found using Eq. (3), which yields a set of equations to be solved simultaneously (see supplementary material for details). In the case of a single component, these equations are
The equations for the twocomponent partition function and Lagrange multipliers can be found in the SI. This method effectively allows us to build a map of the commonality of complex, correlated behaviors relative to basic rates of shape change (as quantified using principal components). For a given Lagrange multiplier governing a particular correlation, a value less than zero indicates a behavior that is less common than expected, and a value greater than zero represents a behavior that is more common.
Stereotypical behavior without biochemical details
After training our model on Dictyostelium shape trajectories, we confirmed that the method had adequately captured the observed correlations by using them to simulate the shape changes of untreated cells responding to cAMP. In order to illustrate the importance of the correlations, we also ran control simulations trained only on the basic rates of increase and decrease in each PC without these correlations. We compared the activities of the uncorrelated and correlated simulations against the observed data. The uncorrelated model acts entirely proportional to the observed rates (though, interestingly, did not match them; Fig. 2D). In contrast, individual cells from the experimental data show very strong anticorrelation, with increases in one component coupled with decreases in the other. This behavior is clearly replicated by the correlated simulations, in both cases appearing in the plot as a red diagonal from the bottom left to the top right. Furthermore, we see suppression of turning behavior in both PCs, with the most poorly represented activity (relative to chance) being a switch in direction in either PC (for example 1+ followed by 1−). This too is reflected in the correlated simulations.
The predictive power of MaxCal simulations goes beyond those correlations on which they were directly trained: We tested the simulations’ ability to predict repetition of any given transition. These patterns took the form of N transitions in T time steps, e.g. five 1+ transitions in ten timesteps. The MaxCal model predicted frequencies of appearance for these patterns that closely resembled the real data (Fig. 3A, model in red, real data in black). In contrast, the uncorrelated model predicted patterns at a much lower rate, for example there are runs of 5 consecutive increases in PC 1 in the real data at a rate of around one in 1.35 minutes. The correlated model predicts this pattern rate to be one in every three minutes. The uncorrelated model predicts the same pattern at a rate of one in 6.67 hours. This result indicates that no higherorder correlations are required to recapitulate the data, allowing us to avoid the huge increase in model complexity that their inclusion would entail.
The greater predictive power of the MaxCal model is reflected by its lower JensenShannon divergence from the observed data for these kinds of pattern (Fig. 3B). The MaxCal model also more closely matches the observed probabilities of generating a given number of transitions in a row, with predictions almost perfect up to 4 transitions in a row (twice the lengthscale of the measured correlations), and far stronger predictive power than the uncorrelated model over even longer timescales (Fig. 3C).
A real world application of MaxCal methods to discriminate between genotypes
We wondered whether the MaxCal methods would accurately discriminate between biologically relevant conditions. To investigate this, we used two comparisons. First, we compared shape data from control AX2 cells against the same cells treated with two drugs targeting cytoskeletal function: the phospholipase A2 inhibitor pbromophenacyl bromide (BPB) and the phosphoinositide 3kinase inhibitor LY294002 (LY) (for details, see^{27}). Second, we compared a stable myosin heavychain knockout against its parent strain (again AX2) (Fig. 4A). We first looked at the effects of these conditions on the distribution of cell shapes, to see whether their effect could be identified from shape, rather than behavior. The drug treatment caused a substantial change in the distribution within the population (Fig. 4B), but still left a substantial overlap. In contrast, the mhcA^{−} cells showed no substantial difference to their parent in shape distribution (Fig. 4C). In both cases, the identification of a condition from a small sample of shape data would not be feasible.
We then compared the behavioral Lagrange multipliers of each condition, found by MaxCal, producing distributions for the estimated values of these by bootstrapping our data (sampling with replacement). The values of \({\gamma }_{\,1\,1}^{+}\) and \({\gamma }_{\,2\,2}^{+}\) are lower in the untreated condition than those in drugtreated condition, indicating the persistence of shape change in WT cells (Fig. 4D). The anticorrelation between PCs 1 and 2 through pseudopod splitting is reflected in \({\gamma }_{\,1\,2}^{+}\) and \({\gamma }_{\,1\,2}^{+}\), both of which have values greater than 1 in WT cells. In comparison, the drugtreated cells have only a moderate anticorrelation. In the mhcA^{−} strain, the differences in the values of \({\gamma }_{\,1\,1}^{+}\), \({\gamma }_{\,1\,2}^{+}\) and \({\gamma }_{\,1\,2}^{+}\) when compared with their parent show similar changes to those observed in drug treatment (Fig. 4E). In both cases, the differences highlighted by these dynamical measurements are striking.
We then applied the MaxCal model to the task of classification. We settled upon classification using knearestneighbors (kNN). In order to see how the strength of our prediction improved with more data, we classified based on the preferred class of N repeats, all known to come from the same data set. We estimated the classification power of our methods by cross validation, dividing the drugtreated data and its control into three sets containing different cells, and dividing the mhcA^{−} and its parent by the day of the experiment. We first performed the classification by shape alone, taking small subsamples of frames from each cell and projecting them into their shape PCs, with our distance measure for the kNN being the Euclidean distance in these PCs. With one, two or 3 PCs, we were able to achieve reasonable classification of the drugtreated cells against their parent as data set size increased, with the accuracy of classification leveling off at around 0.85 (with 1 being perfect and 0.5 being no better than random chance, Fig. 5A–C, blue). In contrast, classification of mhcA^{−} cells was little better than random chance, even with relatively many data (Fig. 5A–C, green). This is unsurprising given the similarity of the distributions of these two conditions. We then calculated our MaxCal multipliers for subsamples of each of these groups, bootstrapping 100 estimates from 20% of each set. We then repeat our kNN classification, instead using for a distance measure the two MaxCal values that best separate our training classes. As the test data come from entirely separated sources (in the case of drugtreated cells coming from different cells, and in the case of the mhcA^{−} being taken on different days), we can be confident that we do not conflate our training and test data. In both the drugtreated case and the mhcA^{−} mutant, the dynamics differ very cleanly between our test and control data. As such, our classification is close to perfect even for only a few samples (Fig. 5D).
As the two Lagrange multipliers that best classified the data both encoded correlations between adjacent timesteps, we guessed that this shortterm memory might be key to recapitulating the dynamic properties of cell shape change. A key aspect of the shape dynamics of AX2 cells is the anticorrelation between the first two PCs at the singlecell level (which is definitionally absent at the population level, as PCA produces uncorrelated axes). To see if memory is vital to recapitualting this dynamical aspect of cell shape change, we constructed two versions of the master equation for our MaxCal system (see SI for details). The first is Markovian (that is, at a given time the probabilities of each possible next event only depend on the current state of the system). We ran Gillespie simulations corresponding to this master equation, and compared the correlations of trajectories from these simulations with those from real data. The expected anticorrelation is clearly observed in the data (Fig. 5E, black line), but the trajectories of our Markovian Gillespie simulations fail to recapitulate it (Fig. 5E, blue line).
We then introduced a memory property to the simulations, allowing the probabilities of each possible event to depend on the nature of the previous event (with the very first event taken using the uncorrelated probabilities). The model has nine possible states (with each state containing its own set of event probabilities), corresponding to the nine possible events that might have preceeded it. These are an increase, a decrease, or no change in each PC independently (3 × 3 events). These nonMarkovian simulations recovered the distribution of correlations observed in the data (Fig. 5E, red line). This indicates that such features of cell shape change can only be addressed by methods that acknowledge a dependency on past events.
Discussion
Eukaryotic chemotaxis emerges from a vast network of interacting molecular species. Here, instead of examining the molecular details of chemotaxis in Dictyostelium discoideum, we have inferred properties that capture cell behavior from observations of shape alone. For this purpose, we quantified shape using Fourier shape descriptors, reduced these shape data to a small, tractable dimensionality by principal component analysis, and built a minimal model of behavior using the maximum caliber variational method. Unlike conventional modeling approaches, such as master equations and their simplifications, our method is intrinsically nonMarkovian, capturing memory effects and shortterm history in the values of the behavioral signature it yields (see SI for further discussion of memory, and a comparison to the master equation). Our approach has the advantage of ease, requiring only the observation of what a cell naturally does^{14}, without tagging or genetic manipulation, as well as of generality, being independent from the specific and poorly understood biochemistry of any one cell type. This is important to understanding chemotaxis, as the biochemistry governing this process can vary greatly: for example, the spermatazoa of C. elegans chemotax and migrate with no actin at all^{28}, but strategies for accurate chemotaxis might be shared among biological systems and cell types.
A number of recent studies have demonstrated the importance of pairwise or shortscale correlations in determining complex behaviors both in space and time. The behavior of large flocks of starlings can be predicted from the interactions of individual birds with their nearest neighbors^{29,30}, and the pairwise correlations between neurons can be used to replicate activity at much higher coupling orders, correctly reproducing the coordinated firing patterns of large groups of cells^{23}. Furthermore, cells in many circumstances use shortrange spatial interactions to organise macroscopically^{31}. Interestingly, these systems appear to exhibit selforganized criticality^{32,33}, in which the nature of their shortrange interactions leads to periods of quiescence puntuated by sudden changes. This could indicate the coupling strengths inherent to a system (such as the temporal correlations in our shape modes in Dictyostelium cells) are crucial for complex behavior. Absence of this behavior could be an indicator of disease as illustrated by both of our aberrant cell types.
Here, we employ a very simple classifier to demonstrate the usefulness of our MaxCal multipliers as a measurement by which we can classify cell behaviours. We choose MaxCal because it is a minimal, statistical approach to modelling a complex phenomenon, allowing high descriptive power with no assumptions made about the underlying mechanism. As our understanding of the molecular biology controlling cell shape improves, an interesting alternative would be to use our data in training recurrent neural network (RNN) autoencoders, a selfsupervised method in which the neural network trains a model to accurately represent the input data. In particular, long shortterm memory RNNs have recently been used to accurately identify mitotic events^{34} in cell video data and classes of random migration based on cell tracks^{35}. The two approaches are not mutually exclusive; MaxCal can provide a neat, compressed basis in which to identify behavioural states of cells, whilst RNNs could be used to learn timeseries rules for transitions between behavioural states.
It is increasingly clear that cell shape is a powerful discriminatory tool^{10}. For example, diffeomorphic methods of shape analysis have the power to discriminate between healthy and diseased nuclei^{36}. Shape characteristics can also be used as an indicator of gene expression^{15}: an automated, shapebased classification of Drosophila haemocytes recently showed that shape characteristics correlate with gene expression levels, specifically that depletion of the tumor suppressor PTEN leads to populations with high numbers of rounded and elongated cells. Of particular note is the observation from this study that genes regulate shape transitions as opposed to the shapes themselves, illustrating the importance of tools to quantify behavior as well as shape. This may be an appropriate approach to take if, for example, creating automated assistance for pathologists when classifying melanocytic lesions (a task which has already proved tractable to automated image analyses^{37}), as classes are few in number, predefined and extensive training data are available. A drawback of the method used by^{15} is that their classes are decided in advance, and the divisions between them are arbitrary. This means that the method cannot find novel important features of shape by definition, as it can only pick between classes decided upon by a person in advance.
A stronger alternative would be to take some more general description of shape and behavior (such as the one we detail here), which could be used to give biopsied cells a quantitative signature. Training would then map these data not onto discrete classes, but onto measured outcomes based on the longterm survival of patients. It will be important for any such method to account for the heterogeneity of primary tissue samples as small subpopulations, lost in gross measurements, may be key determinants of patient outcomes. Such an approach would allow a classifier to identify signs of disease and metastatic potential not previously observed or conceived of by the researchers themselves. As machine learning advances, it will be vital to specify the problem without the insertion of our own biases. Then, behavioral quantification will become a powerful tool for medicine.
Methods
Cell culture
The cells used in our experiments are either of the Dictyostelium discoideum AX2 strain, or a stable myosin heavychain knockout (mhcA^{−}) in an AX2 background. Cells are grown in a medium containing 10 μg/mL Geneticin 418 disulfate salt (G418) (SigmaAldrich) and 10 μg/mL Blasticidine S hydrochloride (SigmaAldrich). Cells are concentrated to c = 5 × 10^{6} cells/mL in shaking culture (150 rpm). Five hours prior to the experiment, cells are washed with 17 mM KNa PBS pH 6.0 (SigmaAldrich). Four hours prior to the experiment, cells are pulsed every 6 minutes with 200 nM cAMP, and are introduced into the microfluidic chamber at c = 2.5 × 10^{5} cells/mL. Measurements are performed with cells starved for 5–7 h. Drugtreated cells were exposed to 200 pM pbromophenacyl bromide and 50 nM LY294002. No. cells sampled in AX2 control, drugtreated, AX2 parent, mhcA^{−} are, respectively, 313, 23, 858, 198.
Microfluidics and imaging
The microfluidic device is made of a μslide 3in1 microfluidic chamber (Ibidi) as described in (12), with three 0.4 × 1.0 mm^{2} inflows that converge under an angle of α = 32° to the main channel of dimension 0.4 × 3.0 × 23.7 mm^{3}. Both side flows are connected to reservoirs, built from two 50 ml syringes (Braun Melsungen AG), separately connected to a customized suction control pressure pump (Nanion). Two micrometer valves (Upchurch Scientific) reduce the flow velocities at the side flows. The central flow is connected to an infusion syringe pump (TSE Systems), which generates a stable flow of 1 ml/h. Measurements were performed with an Axiovert 135 TV microscope (Zeiss), with LD PlanNeofluar objectives 20x/0.50N.A. and 40x/0.75N.A. (Zeiss) in combination with a DV2 DualView system (Photometrics). A solution of 1μM Alexa Fluor 568 hydrazide (Invitrogen) was used to characterize the concentration profile of cAMP (SigmaAldrich) because of their comparable molecular weight.
Image preprocessing
We extracted a binary mask of each cell from the video data using Canny edge detection, thresholding, and binary closing and filling. The centroid of each mask was extracted to track cell movement. Overlapping masks from multiple cells were discarded in order to avoid unwanted contact effects, such as distortions through contact pressure and cellcell adhesion. For each binary mask, the coordinates with respect to the centroid of 64 points around the perimeter were encoded in a complex number, with each shape therefore recorded as a 64 dimensional vector of the form S = x + iy. These vectors were passed through a fast Fourier transform in order to create Fourier shape descriptors. Principal component analysis was performed on the power spectra (with the power spectrum P(f) = s(f)^{2} for the frequencydomain signal s(f)) to find the dominant modes of variation. This approach is superior to simple descriptors such as circularity and elongation, as key directions of variability within the highdimensional shape data cannot be known apriori. As we have previously reported^{13}, 90% of Dictyostelium shape variation can be accounted for using the first three principal components (PCs), corresponding to the degree of cell elongation (PC 1), pseudopod splitting (PC 2) and polarization in the shape (PC 3) (Fig. 1B), with around 85% of variability accounted for in just two, and 80% in one.
References
 1.
Csermely, D., Mainardi, D. & Agostini, N. The predatory behaviour of captive wild kestrel, Falco tinnunculus L. Bull Zool 56, 317–320 (1989).
 2.
Stephens, G. J., JohnsonKerner, B., Bialek, W. & Ryu, W. S. Dimensionality and dynamics in the behavior of C. elegans. PLoS Comput Biol 4, e1000028 (2008).
 3.
Andrew, N. & Insall, R. H. Chemotaxis in shallow gradients is mediated independently of PtdIns 3kinase by biased choices between random protrusions. Nat Cell Biol 9, 193–200 (2007).
 4.
Neilson, M. P. et al. Chemotaxis: A feedbackbased computational model robustly predicts multiple aspects of real cell behaviour. PLoS Biol 9, e1000618 (2011).
 5.
van Haastert, P. J. M. A model for a correlated random walk based on the ordered extension of pseudopodia. PLoS Comput Biol 6, e1000874 (2010).
 6.
Otsuji, M., Terashima, Y., Ishihara, S., Kuroda, S. & Matsushima, K. A conceptual molecular network for chemotactic behaviors characterized by feedback of molecules cycling between the membrane and the cytosol. Sci Signal 152, ra89 (2010).
 7.
Westendorf, C. et al. Actin cytoskeleton of chemotactic amoebae operates close to the onset of oscillations. Proc Nat Acad Sci USA 110, 3853–3858 (2013).
 8.
Davidson, A. J, Amato, C., Thomason, P. A. & Insall, R. H. WASP family proteins and formins compete in pseudopod and blebbased migration. J Cell Biol jcb.201705160 (2017).
 9.
Swaney, K. F., Huang, C. H. & Devreotes, P. N. Eukaryotic chemotaxis: a network of signaling pathways controls motility, directional sensing, and polarity. Ann Rev Biophys 39, 265–289 (2010).
 10.
Marklein, R. A., Lam, J., Guvendiren, M., Sung, K. E. & Bauer, S. R. FunctionallyRelevant Morphological Profiling: A Tool to Assess Cellular Heterogeneity. Trends Biotech 36, 105–118 (2017).
 11.
van Haastert, P. J. M. & Postma, M. Biased random walk by stochastic fluctuations of chemoattractantreceptor interactions at the lower limit of detection. Biophys J 93, 1787–1796 (2007).
 12.
Tweedy, L., Knecht, D. A., Mackay, G. M. & Insall, R. H. SelfGenerated Chemoattractant Gradients: Attractant Depletion Extends the Range and Robustness of Chemotaxis. PLoS Biol 14(3), e1002404 (2016).
 13.
Tweedy, L., Meier, B., Stephan, J., Heinrich, D. & Endres, R. G. Distinct cell shapes determine accurate chemotaxis. Sci. Rep. 3, https://doi.org/10.1038/srep02606 (2013).
 14.
Keren, K. et al. Mechanism of shape determination in motile cells. it Nature 453, 475–480 (2008).
 15.
Yin, Z. et al. A screen for morphological complexity identifies regulators of switchlike transitions between discrete cell shapes. Nat Cell Biol 15, 860–872. (2013).
 16.
Broekmans, O. D., Rodgers, J. B., Ryu, W. S. & Stephens, G. J. Resolving coiled shapes reveals new reorientation behaviors in C. elegans. eLife 5, e17227 (2016).
 17.
Gyenes, B. & Brown, A. E. X. Deriving ShapeBased Features for C. elegans Locomotion Using Dimensionality Reduction Methods. Front Behav Neurosci 10, 159 (2016).
 18.
Valletta, J. J., Torney, C., Kings, M., Thornton, A. & Madden, J. Applications of machine learning in animal behaviour studies. Animal Behav. 124, 203–220 (2017).
 19.
GomezMarin, A., Stephens, G. J. & Brown, A. E. X. Hierarchical compression of Caenorhabditis elegans locomotion reveals phenotypic differences in the organization of behaviour. Front Behav Neurosci 13, 20160466 (2016).
 20.
Pressé, S., Ghosh, K., Phillips, R. & Dill, K. A. Dynamical fluctuations in biochemical reactions and cycles. Phys Rev E 82, 031905 (2010).
 21.
Pressé, S., Ghosh, K., Lee, J. & Dill, K. A. Principles of maximum entropy and maximum caliber in statistical physics. Rev Mod Phys 85, 1115–1141 (2013).
 22.
Jaynes, E. T. The minimum entropy production principle. Ann Rev Phys Chem 31, 579–601 (1980).
 23.
Schneidman, E., Berry, M. J., Segev, R. & Bialek, W. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440, 1007–1012 (2006).
 24.
Cavagna, A. et al. Dynamical maximum entropy approach to flocking. Phys Rev E 89, 042707 (2014).
 25.
Vasquez, J. C., Marre, O., Palacios, A. G., Berry, M. J. II & Cessac, B. Gibbs distribution analysis of temporal correlations structure in retina ganglion cells. J PhysiologyParis 106, 120–127 (2012).
 26.
Firman, T., Balázsi, G. & Ghosh, K. Building Predictive Models of Genetic Circuits Using the Principle of Maximum Caliber. J Biophys J 113, 2121–2130 (2017).
 27.
Meier, B. et al. Chemotactic cell trapping in controlled alternating gradient fields. Proc Natl Acad Sci USA 108, 11417–11422 (2011).
 28.
Nelson, G. A., Roberts, T. M. & Ward, S. Caenorhabditis elegans spermatozoan locomotion: amoeboid movement with almost no actin. J Cell Biol 92, 121–131 (1982).
 29.
Toner, J. & Tu, Y. Longrange order in a twodimensional dynamical XY model: how birds fly together. Phys Rev Lett 75, 4326–4329 (1995).
 30.
Bialek, W. et al. Statistical mechanics for natural flocks of birds. Proc Nat Acad Sci USA 109, 4786–4791 (2012).
 31.
De Palo, G., Yi, D. & Endres, R. G. A criticallike collective state leads to longrange cell communication in Dictyostelium discoideum aggregation. PLoS Biol 15, 1–25 (2017).
 32.
Mora, T. & Bialek, W. Are biological systems poised at criticality? J Stat Phys 2, 268–302 (2011).
 33.
Chialvo, D. R. Emergent complex neural dynamics. Nat Phys 6, 774–750. (2010).
 34.
Phan, H. T. H., Kumar, A., Feng, D., Fulham, M. & Kim, J. Unsupervised twopath neural network for cell event detection and classification using spatiotemporal patterns. IEEE T Med Imaging, https://doi.org/10.1109/TMI.2018.2885572 (2018).
 35.
Kimmel, J. C, Brack, A. S. & Marshall, W. F. Deep convolutional and recurrent neural networks for cell motility discrimination. bioRxiv, https://doi.org/10.1101/159202 (2019).
 36.
Rohde, G. K., Ribeiro, A. J. S., Dahl, K. N. & Murphy, R. F. Deformationbased nuclear morphometry: capturing nuclear shape variation in HeLa cells. Cytom Part A 73A, 341–350 (2008).
 37.
Esteva, A. et al. Dermatologistlevel classification of skin cancer with deep neural networks. Nature 542, 115–118 (2017).
Acknowledgements
We are grateful to Börn Meier for sharing his data, and to both André Brown, Linus Schumacher and Peter Thomason for a critial reading of the manuscript. This work was supported by Cancer Research UK core funding (L.T.), the Deutsche Forschungsgemeinschaft (DFG fund HE5958/21), the Volkswagen Foundation grant I/85100 (D.H), and the BBSRC grant BB/N00065X/1 (R.G.E.) well as the ERC Starting Grant 280492PPHPI (R.G.E.).
Author information
Affiliations
Contributions
L.T. and R.G.E. designed the study. L.T. and P.W. performed the experiments, and L.T. conducted data analysis and modelling. All authors (L.T., P.W., D.H., R.H.I., R.G.E.) analyzed results and data, and wrote the paper.
Corresponding author
Ethics declarations
Competing Interests
The authors declare no competing interests.
Additional information
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Tweedy, L., Witzel, P., Heinrich, D. et al. Screening by changes in stereotypical behavior during cell motility. Sci Rep 9, 8784 (2019). https://doi.org/10.1038/s4159801945305w
Received:
Accepted:
Published:
Further reading

Physicsinformed deep learning characterizes morphodynamics of Asian soybean rust disease
Nature Communications (2021)

Artificial Intelligence in Intracoronary Imaging
Current Cardiology Reports (2020)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.