Analysis of protrusion dynamics in amoeboid cell motility by means of regularized contour flows

Amoeboid cell motility is essential for a wide range of biological processes including wound healing, embryonic morphogenesis, and cancer metastasis. It relies on complex dynamical patterns of cell shape changes that pose long-standing challenges to mathematical modeling and raise a need for automated and reproducible approaches to extract quantitative morphological features from image sequences. Here, we introduce a theoretical framework and a computational method for obtaining smooth representations of the spatiotemporal contour dynamics from stacks of segmented microscopy images. Based on a Gaussian process regression we propose a one-parameter family of regularized contour flows that allows us to continuously track reference points (virtual markers) between successive cell contours. We use this approach to define a coordinate system on the moving cell boundary and to represent different local geometric quantities in this frame of reference. In particular, we introduce the local marker dispersion as a measure to identify localized membrane expansions and provide a fully automated way to extract the properties of such expansions, including their area and growth time. The methods are available as an open-source software package called AmoePy, a Python-based toolbox for analyzing amoeboid cell motility (based on time-lapse microscopy data), including a graphical user interface and detailed documentation. Due to the mathematical rigor of our framework, we envision it to be of use for the development of novel cell motility models. We mainly use experimental data of the social amoeba Dictyostelium discoideum to illustrate and validate our approach.


Introduction
Amoeboid motion is one of the most widespread forms of cell motility in the living world [1].It plays a key role in many essential functions of the human body, such as responses of the immune system [2] or the healing of injured tissue [3].Its medical relevance also extends to the field of cancer research, as metastatic tumor cells rely on amoeboid motility to invade the surrounding tissue [4].Amoeboid locomotion is based on dynamical changes of the cell shape.Specifically, localized protrusions of the cell membrane, often called pseudopodia, that are extended in the direction of motion, are generally seen as the basic morphological entities that drive amoeboid movement [5].Together with membrane contraction at the back of the cell body, their extension results in a displacement of the center of mass of the cell.This requires, besides the coordinated pattern of protrusion and retraction, also the formation and rupture of adhesive contacts to a substrate or to a surrounding extracellular matrix [6].
The mechanical forces that drive the shape changes of amoeboid cells are generated by the actin cytoskeleton, a dynamic filament meshwork at the inner face of the cell membrane [7].The main building blocks of this network are actin filaments that are subject to a constant turnover by polymerization and depolymerization, resulting in a continuous rapid reorganization of the network structure.This process is assisted by a host of auxiliary cytoskeletal proteins that initiate the nucleation, capping, or severing of actin filaments as well as filament bundling, branching, and membrane cross-linking.The cytoskeletal machinery is orchestrated by biochemical signaling pathways that coordinate the spatiotemporal patterns of activity in the actin system across the cell cortex [8].These upstream signaling pathways also provide a link to membrane receptor signals, so that cells may react to gradients of extracellular cues by moving directionally towards a chemical source-a process commonly referred to as chemotaxis [9].
The mathematical modeling of amoeboid motion is a long-standing challenge that has been addressed at different levels of complexity, ranging from random walk models for the center of mass of the cell [10][11][12] up to detailed high-dimensional models for the intracellular signaling activity [13].Typically, current models focus on selected mechanistic aspects of amoeboid motility and describe, for example, actin dynamics, cell-to-cell variability, or the switching between different migratory modes in more detail, see for example [14][15][16].Also, first attempts have been made to incorporate several key components, such as dynamic signaling patterns, polarity formation, and cytoskeletal activity, in a modular approach [17].These models, however, remain qualitative and their comparison with experimental data oftentimes relies on visual inspection.As the entire biological process involves many hundreds of interacting proteins and signaling molecules, with many mechanistic details yet unknown, a quantitative model that includes the full molecular details remains out of reach.
To advance the mathematical description of amoeboid motility, we envision that current efforts of mechanism-based modeling are complemented by a more systematic, data-driven approach.This requires a mathematical framework that allows us to systematically develop a quantitative model based on experimental data.Such a framework should rely on observables that encode the key characteristics of amoeboid motility and, at the same time, are readily accessible experimentally.Trajectories of the center-of-mass of the cell can be easily recorded in large amounts from September 9, 2021 2/31 low-resolution bright-field microscopy data but reflect only very limited, integral information on the entire process.The intracellular signaling pathways and the cytoskeletal mechanisms, on the other hand, are difficult to access and knowledge on this part of the system remains highly incomplete.We therefore concentrate on the cell shape as the central reference quantity.The cell shape is fully accessible by standard microscopy techniques and can be easily recorded with sufficient spatial and temporal resolution.Moreover, its dynamic evolution implicitly reflects the intracellular processes and determines the center-of-mass trajectory of the cell.
In our long-term quest for a quantitative, data-driven model of amoeboid motility, several steps are required: First, the development of a mathematical framework for the description of experimentally observed shape dynamics; second, the design of a model of the contour dynamics that predicts realistic shape evolutions; and finally, the incorporation of mechanistic information on key intracellular processes as the driving determinants of the contour dynamics.These become accessible by imaging of fluorescently tagged fusion proteins and by more advanced methods, such as knock-sideways and optogenetic approaches.Here, existing mechanistic models may provide a useful basis and might merge into a joint modeling concept.
In this article, we concentrate on the first aspect of the above agenda.We provide a mathematically well-defined approach that allows for a detailed analysis of the complex, multifaceted contour dynamics of amoeboid cells.A key ingredient is the concept of regularized flows between contours that define an evolution of virtual markers in time.Using contour flows, we define a coordinate system on the evolving contours (strongly regularized case) and approximate local quantities of interest (weakly regularized case).While the strongly regularized flow is used to define trajectories of virtual markers over the entire time span, the weakly regularized flow is used to obtain information on local membrane changes, which is mapped subsequently onto the global flow representing the coordinate system.This separation between global and local flows is an essential feature of our approach.
Established approaches to describe virtual markers on an evolving contour most commonly include level set methods (LSM), e.g., in [18], where the cell contour dynamics was additionally decomposed into a translation of the entire cell and the deformation of its contours.In [19], the LSM was compared to a mechanistic spring model penalizing a dense concentration of virtual markers.While the mechanistic model provides better computation times than the LSM, it shows limitations regarding topological mapping violations during strong shape deformations.Alternatively, electrostatic field equations were successfully employed, tackling mapping violation issues while providing better computation times [20].However, it does not provide trajectories of virtual markers over the entire time interval in its current form.In [21], a mapping was chosen which minimizes the sum of squared distances between virtual markers (so-called mean squared displacement), while enforcing an evenly-spaced distribution of virtual markers.Our work aims at (i) presenting a theoretical framework for the analysis of marker dynamics and (ii) combining the respective merits of the above approaches into a single method.
In the spirit of previous approaches, we rely on the widely used concept of kymographs to graphically represent the space-time dynamics of different geometric quantities, such as the speed of membrane displacement or the local curvature, along the cell contour.In this context, we propose a novel criterion based on a single dynamic quantity for defining membrane expansions/contractions.Prior approaches to identify membrane protrusions relied on the simultaneous matching of multiple criteria, such as critical values of curvature, protrusion speed, and pseudopod lifetime [22] or identified protrusion events as points in time only [21], without providing major protrusion properties, such as area growth rate, shape, and others.
In [23], pseudopods were detected by using a hierarchical clustering algorithm in which individual membrane extensions are connected based on direction and continuity in space and time.Furthermore, image skeletonization on contours was used to identify and characterize pseudopods in an automated way [24].Unfortunately, all existing approaches have in common an undesirable blending of the protrusion-defining criteria with their numerical implementation.This makes it difficult to discern biological effects from numerical artifacts.We instead first define our criterion in mathematical terms and only subsequently implement it numerically, allowing us to control numerical errors and avoid artifacts.
The methodology of this article is implemented in AmoePy, a Python-based toolbox for analyzing and simulating amoeboid cell motility [25].Furthermore, AmoePy features multiple data handling and analysis tools, an easy-to-use graphical user interface, and data files containing multiple experimental cell tracks in addition to simple artificial tracks which were used to additionally validate the algorithm.
The article is organized as follows.We first develop the mathematical framework and introduce the concept of regularized contour flows.We then illustrate our approach in applications to experimental recordings of the social amoeba Dictyostelium discoideum, a widely used model organism for the study of amoeboid motility.Finally, based on the results of our analysis, we illustrate in the Discussion section how we envision the next steps towards a quantitative, data-driven model of amoeboid motility, based on a point process of the protrusive activity.

Data acquisition and image segmentation
Experiments were performed with AX2 cells of the social amoeba Dictyostelium discoideum.As a marker for filamentous actin, cells expressed fluorescently tagged Lifeact (C-terminally fused with mRFP, plasmid kindly provided by Igor Weber, Zagreb, HR).They were grown at 20°C in liquid culture flasks containing HL5 medium including glucose (Formedium, Hunstanton, GB) and 10 µg/ml G-418 sulfate (Cayman Chemical Company, US) as a selection marker.Before each experiment, cells were harvested from culture flasks and grown in 25 ml overnight shaking cultures at 180 rpm under otherwise identical conditions.Afterward, nutrients were removed by centrifugation and washing of the cell pellet with Sørensen phosphate buffer at pH 6 (14.7 mM KH 2 PO 4 , Merck, Darmstadt, DE; 2 mM Na 2 HPO 4 × H 2 O, Merck, Darmstadt, DE).Then, cells were resuspended in fresh buffer and droplets were formed in a Petri dish to initiate the streaming processes.
For image acquisition cells were transferred after 5 hours to a glass bottom dish (Fluorodish, ibidi GmbH, Martinsried, DE).During imaging, they were kept in Sørensen phosphate buffer at 20°C.Images were taken with a Zeiss LSM780 laser scanning confocal microscope (Carl Zeiss AG, Oberkochen, DE) at a frame rate of one image per second, using a 63× or 40× oil immersion objective.Fluorophores were excited at 651 nm and emission was recorded between 562 nm and 704 nm.For details see also [15], where the same data set was used in a different context.
Fluorescence image (8-bit gray scale) were segmented using a modified version of the active contour (snake) algorithm described in [21,26].Based on this algorithm, we parameterized the cell boundary in each frame by a closed string of M = 400 equidistant nodes.As frames were taken at discrete time points t 0 , . . ., t K−1 with equal time difference δt = t k+1 − t k = 1 sec, we denoted the discrete representation of the cell contour in a given frame at time t k as γ k with supporting points (1) September 9, 2021 4/31 In the result section, the data sets consists of K = 500 to 1000 time frames.For later reference, we set t 0 = 0 and t K−1 = T .

Estimate of contour dynamics
We used a real-valued smoothing spline for the x and y coordinates based on Gaussian process regression (GPR) using a Poisson kernel, for details see S1 Text.This yielded a parametrization of the contour Γ k at time frame in terms of a finite sum of smooth kernels (see e.g.[27,28] for details) where P is a suitably scaled Poisson kernel.Support points were chosen to correspond to normalized secant length along the contour for m = 0, . . ., M − 1 and k = 0, . . ., K − 1, see S1 Text for details.We parametrized the contour in the mathematical positive sense, i.e., the interior of the cell is on the left when going around the contour with increasing θ.In the numerical implementation, we used the rescaled arc length coordinates, which we denoted again by θ.This gives denoting the length of contour Γ k .Note that Φ k is only uniquely determined up to a phase shift, i.e., for every Φ k and τ , also Φ k,τ (•) = Φ k (• + τ ) is a valid parametrization of Γ k .The phase shift was chosen by an additional requirement in the next section.The smooth parametrization Φ k allowed us compute local quantities along the contour, e.g., its curvature where R π/2 is an anti-clockwise rotation by π/2.As global quantity, we determined the center of mass C k of contour Φ k as Since the segmentation points were rather densely spaced over the contour, they well constrained the smooth contours.Based on the kernel representation, all geometric quantities, such as arc length and curvature were defined intrinsically for each contour, and may be easily computed numerically and with high precision.Connecting the contours along the time axis, however, is not intrinsically well defined, and is bound to choices.In a first step, we constructed a global coordinate flow, which served as a reference frame for further local flows to be defined subsequently.In the limit of infinitely densely sampled contours, this global coordinate system results in a mapping Vice versa, any coordinate system defines a coordinate flow over the tube of contours, i.e., the cell contours in 2d mapped into a 3d space-time coordinate system.If p 0 = (x 0 , y 0 ) is a point on the first contour at t = 0 and if θ 0 is its arc length coordinate, then t → Φ(t, θ 0 ) corresponds the movement of p 0 over the space-time tube of contours, see S1 Fig. Additionally, in panel (B), the velocity v(θ) that propagates a virtual marker θ over a particular distance to the next contour, is illustrated as well as its decomposition into an outward-pointing normal velocity v n (θ) and a tangential velocity v t (θ).
For the weakly and non-regularized cases, a stretching effect of virtual markers can be observed for expanding parts (gray area), whereas clustering effects of virtual markers occur at contracting regions.Assuming infinitesimally small time steps, the "stretching rate" at a virtual marker θ that arises from transitioning from one contour September 9, 2021 6/31 to the next one is given by where v t and v n denote the (scalar) speed in tangential and normal direction, respectively, and κ the curvature.In the following sections, we show that this kind of "stretching rate", in the sequel termed instantaneous dilation rate, is a useful quantity to identify active regions of the contour during cell motility.
While the normal component of the velocity can be approximated easily by taking the normal distances from the first contour to the second divided by δt, the tangential component remains unknown and therefore requires additional attention.Naive flows, such as the shortest path flow from Fig 1 (right column) or the flow merely based on the normal component, can produce topological mapping violations.An example of such a mapping violation can be seen in panel (F) where the order of virtual markers on the second contour is inconsistent with the previous one.
In general, additional constraints on the tangent component are necessary with the aim of preserving an evenly-spaced distribution of virtual markers over time.In [19], these additional constraints were formulated as a mechanistic spring model, providing better computation times than the more commonly used LSM.However, it is mentioned that for strong shape deformations this mechanistic model easily fails because of mapping violations.This approach, as well as ours, is based on a cost functional minimizing the trajectories of virtual markers to the next contour while penalizing the distances of neighboring virtual markers.While in the mechanistic spring model the penalizing distance is measured in R 2 , the space of contours, our approach measures the penalizing distance in [0, 2π), the space to parametrize the contour.This key difference resolves the issue of mapping violations during large deformations in our approach.
In S2 Fig, contour mappings for each of these two regularization schemes are displayed showing the advantageous stability of our developed S 1 regularization during membrane spikes and other large shape deformations.In S3 Fig, a selection of contour mappings between two frames with larger shape deformations are displayed.The underlying microscopy data was recorded with an imaging rate of δt ≈ 3.13s.

Our method combines and improves different aspects of previous approaches:
• Spatio-temporal tracking of a fixed number of virtual markers, • Distances between virtual markers flexibly adjustable in terms of a single regularization, • Faster than computationally expensive approaches such as LSM, • Capability to capture large shape deformations/avoidance of mapping violations, • Simplicity regarding the interpretation & number of parameters.
Another key feature of our approach is the usage of two different flows, separating the underlying coordinate system, defined by one (strongly regularized) flow, and the dynamic quantities of interest, which are defined separately, e.g., by a weakly regularized flow.
The maximal correlation coordinate system (MCCS).The starting point are the parameterized contours Φ k in Eq (2) for k = 0, . . ., K − 1.To make the influence of the sampling rate more prominent, we also used the notation September 9, 2021 7/31 As stated above, the parametrization of Φ k is only determined up to a phase shift by Eq (3).We finally chose the phase shift and therefore the parametrization of Γ k+1 by minimizing the distance to the previous contour Γ k in a mean squared sense, i.e., We may alternatively interpret Eq (7) as optimizing the cross covariance between the two contours when interpreted as vector-valued functions In the sequel, we used Φ k+1 (•) = Φ k+1 (• + τ k+1 ) and omitted the tilde for ease of notation.Effectively, choosing the phase shift amounts to fixing the 'zero point' Φ k+1 (0) on Γ k+1 .This is the coordinate system that from now on is used to represent the contour geometry, i.e., each scalar quantity q defined on the contour Γ k with q = q(t k , (x, y)) and (x, y) ∈ Γ k as function (k, θ) → q(kδt, θ) of discrete time and continuous space can be represented w.r.t. the chosen coordinate system Φ.
The Eulerian and Lagrangian points of view.Any flow that maps the contour Γ k into Γ k+1 is determined by a mapping which describes the translation along the contour φ k .To ensure that θ → φ k (θ) is a one-to-one map, we required in addition An example of a mapping violation, i.e. a position θ with , where the order of virtual markers on the following contour is inconsistent with the previous one.The iteration describes the trajectory (θ k ) k=0,...,K−1 of the starting point at coordinate θ 0 on the first contour in our coordinate system.This approach to visualize the flow shall be called the Eulerian point of view, since it describes the translation vector field from Γ k to Γ k+1 in the coordinate system of Γ k : The Lagrangian point of view instead describes the flow in the coordinate system it generates, which is different from our MCCS.Denote the coordinate of a point on the first contour by its angle coordinate θ 0 , and let χ k be the mapping of Γ 0 to Γ k recursively defined by Both points of view are useful to understand and describe a flow over the contour.The translation vector W k in Lagrangian coordinates is simply and is linked to the Eulerian description via We used the functions φ k , V k and W k interchangeably to specify the flow from Γ k to Γ k+1 .
September 9, 2021 8/31 Transport along the flow.For any flow, we may define the instantaneous dilation rate LD of the flow.Considering two infinitesimally nearby points on contour Γ k , we see that the local relative dilation/contraction factor is obtained from Note that our global coordinate system MCCS induces a flow with uniform dilation rate.To describe the transport of a density under the flow, consider points on the contour Γ k that are distributed according to a density µ k (θ)dθ.Under the flow induced by φ k this density changes according to Starting from µ 0 (•) ≡ 1/(2π), this defines the transported density on all contours.In the Lagrangian picture this transport preserves the density µ 0 (•), which follows from the fact that by definition the starting angle does not change under the flow.The density µ k can be written in Lagrangian coordinates as A regularizing family of flows.We next defined a family of local mappings φ k between successive contours that yields a compromise between the shortest path flow and the uniform dilation coordinate flow, presented in the two end-member cases in Fig 1 .Another suitable name for shortest path flow is reversed normal flow since the incoming trajectories under this flow are always orthogonal to the successive contour.
The mean squared velocity of the flow (with respect to a density µ k ) is given as The normal flow from contour Γ k to Γ k+1 is the flow that departs from the first contour in the normal direction until it intersects with the second contour.The normal flow from Γ k+1 to Γ k shall be called the reversed normal flow from Γ k to Γ k+1 .This is the unconstrained minimizer of F k .If there are no intersections of flow lines, it defines a one-to-one mapping between the two contours.In general, however, direct minimization of the functional F k does not yield a valid flow because of self-intersections, and as a consequence, the induced map is multiple-valued.We therefore need to regularize the flow.A natural requirement is that the flow tends to enforce non-uniformly distributed points on contour Γ k towards more uniformly distributed points on contour Γ k+1 .We proposed to quantify the degree of non-uniformity of a distribution µ(θ)dθ by means of Since any distribution satisfies 2π 0 µ(θ) = 1 and µ(θ) > 0, the minimizer actually corresponds to the uniform distribution.Other measures of non-uniformity are possible, for instance the entropy This kind of regularization has been proposed in [29] by Otto, where the optimal flow is understood as a gradient flow of the entropy potential with respect to the September 9, 2021 9/31 Wasserstein transport distance of the marker density.Also very appealing, in this paper, we used the characterization in Eq (14), since it leads to a more readily tractable numeric quadratic minimization problem.In terms of the defining mapping φ k , the functional U reads as follows from Eq (11).The regularized flow is defined as the flow that minimizes a compromise between both cost functions Note that ] depends on both, φ k and the measure µ k .When iterating over all contours, one needs to update the measure before optimizing the flow for the next time step.There are two end-member cases for the regularized flow: • For large λ the optimal flow essentially immediately uniformizes the density of the initial contour.Thereafter, it is the uniform stretching flow that minimizes the mean square distances between the contours.This is precisely the coordinate flow defined before.
• For small λ the optimal flow allows for arbitrary local stretching rates to minimize the point-wise distance.Here the limit defines the regularized shortest path flow, respectively, the regularized reverse normal flow.
Note that straightforward pointwise minimization of the flow distance from Γ k to Γ k+1 may lead to overlapping connections and hence singular mappings between the contours.If instead, we regularize with small λ, such overlaps are avoided.
The virtual marker picture.For numerical implementation, we discretized the cost functionals using the concept of virtual markers on the contours.The virtual markers are a discretized version of the Lagrangian coordinates.Since µ k is the transported density, the first cost functional for φ k in the Lagrangian point of view using Eq (12) is given by while the second functional using Eq (12) is given by See S1 Text for details of the derivation.Both equations are well suited for a discrete numerical approximation for a given function f and initially equidistant θ 0,i = 2πi/N with i = 0, . . ., N − 1 using If we now approximate the continuous mapping φ k by its discrete values on the virtual marker points θ k = (θ k,0 , . . ., θ k,N −1 ) with September 9, 2021 10/31 then the first cost function may be approximated as and the second cost function as For the entropy based measure of uniformity, consider a collection of points θ k,0 , . . ., θ k,N −1 ∼ µ k on Γ k that are distributed according to the density µ k (•) (not necessarily uniform).Then for any function f , it is In this discrete virtual marker approximation, the local dilation rate at θ i,k , also called the local dispersion, reads Please note that for infinitesimally small time steps another formulation of this rate is given by Eq (6).Especially for level set methods, this formula is more applicable.Finally, the discrete optimization problem is given by Note that in the discrete optimization problem we do not pose any requirements on the space of transformations φ k ensuring condition Eq (8).
Marker re-initialization for weakly-regularized contour flows.In general, the distribution of virtual markers θ k ∼ µ k on Γ k depends on the initial distribution of virtual markers θ 0 ∼ µ 0 on Γ 0 , since the density µ k results from the transport of µ 0 by the flow.As a consequence, functions derived from the local contour flow like, e.g., the local dilation rate along Γ k , may depend on the initial distribution on Γ 0 .By re-initializing the distribution of virtual markers on Γ k for any k = 1, . . ., K − 1, this dependence may be removed.A natural choice is to re-initialize with the uniform distribution, i.e., using with ξ i = 2πi for i = 0, . . ., N − 1 instead of Eq (18).To approximate the local flow φ = (φ k ) k=0,...,K−1 between contours, we thus solved the re-initialized optimization problem φ k,λ = argmin with ξ = (ξ i ) i=0,...,N −1 .For the local dispersion, e.g., this resulted in with no dependence on the initial distribution of markers on Γ 0 .
2. Determine contour based quantities like the curvature in Eq (4) or the center of mass in Eq (5).
3. To determine the (strongly-regularized) coordinate flow, consider N equally spaced markers N on the initial contour and choose a large λ value.Iteratively solve the regularization problem in Eq (24) to determine the coordinate markers θ k+1 for Γ k+1 from the coordinate markers θ k for Γ k .Since both θ k+1 and θ k are approximately equally spaces, solving the minimization problem amount to choosing θ k+1,0 .4. To determine the (weakly regularized) local flow φ = (φ k ) between successive contours, choose a small λ value and solve the regularization problems in Eq (24) to determine the coordinates θ k+1,i = φ k,λ (ξ i ) on Γ k+1 based on N equally spaced markers ξ 0 , . . ., ξ N −1 on Γ k .
5. Determine contour flow based quantities like the local dispersion in Eq (25) or the local motion: see also Eq (9).
6. Map all quantities, based on contour flow as well as contour only, onto the (strongly regularized) global flow.
Numerical implementation and reproducibility.All methods presented in this article are fully accessible as an open source Python-based toolbox AmoePy [25].Implementations of the above methods and additional routines necessary to reproduce the figures and results of this article are part of this toolbox.Furthermore, AmoePy contains: • An object-oriented analysis tool to handle cell contours, i.e., to shorten, extend, or manipulate existing data sets and to extract additional geometric quantities (e.g.area, perimeter, and normal vectors along the cell contour), • A python class to perform Gaussian process regressions in 2d (space) and 3d (space-time) with different selectable kernel functions, • Several routines generating videos of cell tracks with corresponding kymographs and occurring expansion/contraction patterns, • A detailed documentation generated by the Python documentation tool Sphinx, • A graphical user interface able to compute and present kymographs September 9, 2021 12/31 Algorithm 1: RegFlow Input: γ 0 , . . ., γ K−1 , r ∈ (0, 1) , σ n > 0, N ∈ N, λ ≥ 0 Output: Γ 0 , . . ., Γ K−1 , θ eval 0 , . . ., θ eval K−1 /* Initialize evenly distributed starting points */ θ init = NormalizedArcLengthCoordinates(γ 0 ); Algorithm to compute regularized flows.The algorithm input consists of three parameters r, σ n , and N regarding the Gaussian process regression, a regularization parameter λ and the initial (segmented) contours γ 0 , . . ., γ K−1 .The regularized flow is described by the output variables Γ 0 , . . ., Γ K−1 denoting smooth contours evaluated at a finite number of coordinate markers θ eval 0 , . . ., θ eval K−1 .
AmoePy, as well as its graphical user interface, is updated regularly.Future outcomes of our ongoing research regarding for example cell segmentation routines and a forward model to simulate amoeboid cell motility will be also added to AmoePy.
The algorithm starts with initializing equally spaced markers, representing the starting points of the flow (see also the algorithmic workflow).Subsequently, the optimization problem in Eq (24) is solved contour-wise by gradient descent.Here, it is beneficial that the GPR also provides these gradients as a by-product, which decreases In contrast to level set and electrostatic methods, the optimization approach presented here is not relying on the computation of intermediate field lines or contours for sufficiently small grid sizes.Notably, level set methods are computationally more expensive because of the time integration of virtual markers along these field lines [19,20].On the other hand, our algorithm contains additional steps such as the re-initialization of virtual markers via Gaussian process regression and the computation of kymograph quantities which may lead to slower computation times than other empirical mapping algorithms.A direct comparison of computation times of these methods is rather difficult since some algorithms rely on a varying number of virtual markers or a more image-based implementation (e.g. in ImageJ).Moreover, the computation time depends on varying configurations of these methods, e.g., the resolution of the underlying grid in LSM, or the number of iterations/tolerance parameters used in optimization approaches such as ours.Unfortunately, a "biological true" mapping does not exist, with which the accuracy of each of those methods can September 9, 2021 13/31 Algorithm 2: ComputeKymographs  For a cell track of 500 contours and 400 virtual markers, the computation time for generating the mapping between two successive contours was measured for different regularization parameters λ.Not surprisingly, the computation time decreases for the end-member cases λ → 0 and λ ≫ 0, where the underlying functional is only defined by one leading term.For each choice of λ, the median of the computation time is below 2s where most of the cases required sub-second computation times.Furthermore, no mapping violations occurred during the cell track despite larger shape deformations and membrane spikes.Therefore, our method has the potential to be used specifically for cell motility models that rely on a fast computation of stable virtual marker trajectories.

Degree of regularisation controls distribution of virtual markers
We applied our approach to time-lapse microscopy data of the social amoebae D. discoideum.are shown.In contrast to Fig 4, the underlying snake of segmented points was obtained from noisy fluorescence image data resulting in many abrupt changes of the contour's direction.In this context, the curvature for different contour estimates is displayed, highlighting the effect of underfitting and overfitting.Notably, the GPR provides an automated way of estimating the hyperparameters balancing the model's complexity and its error residuals, see S1 Text for more information.This way, an accurate approximation of the cell contour can be obtained even for noisy data while preserving the main characteristics of its shape.
In S8 Fig, the effect of different imaging frequencies on the resulting kymographs is shown.In this context, local motion kymographs are computed by using (i) the entire microscopy data (one image/contour per second) and (ii) down-sampled data sets based on every 2nd, 3rd, 5th, and 10th image/contour.While global characteristics of the cell track are captured even for a lower temporal resolution (δt > 3), the identification of local membrane changes becomes impossible.We therefore recorded each cell track with an imaging rate of one frame per second (δt = 1).See also S3 Fig, in which contour mappings for larger shape deformations are displayed.Here, the underlying cell track was recorded with an imaging rate of δt ≈ 3.13s.
The continuous representations Φ 0 , . . ., Φ K−1 of all contours are the input to the optimization problem Eq (16) to determine the regularized contour flow.for two illustrative contours).In the absence of any regularization (panel (A), λ = 0), virtual markers are thinned out in some regions (linked to expanding areas), while they are clustered in others (in particular at the back of the cell).The case λ = 0 corresponds to minimizing the translation from the first contour to the next, i.e., each virtual marker on the first contour is linked to its nearest neighbor on the second contour.This may result in mapping violations θ k+1,i+1 − θ k+1,i ≤ 0 on the second contour.An example of mapping violations can be seen on the left-hand side in panel (A).As mentioned earlier, the flow obtained with λ = 0 is defined by the shortest path flow or equally the reversed normal flow.Thus, instead of taking the trajectories to the nearest neighbors on the next contour, one could also choose the shortest normal vectors from the second contour to the first one.
In the weakly regularized case, virtual marker thinning and clustering is still prominent (see panel (B) with λ = 0.1), but to a lesser extent.Moreover, in the presence of regularization, virtual marker trajectories are interdependent, which results in trajectories without mapping violations.In the limit of strong regularization, the marker points remain uniformly distributed on every contour, while minimizing the overall distance between contours (see panel (C) with λ = 1000).This makes the strongly regularized virtual marker trajectories an ideal candidate for a time-evolving reference frame and corresponds to the previously defined MCCS coordinate system.
By choosing a strong regularization for the coordinate system, however, September 9, 2021 16/31 information on local contour changes is largely lost.Therefore, we used the strongly regularized case only to determine the coordinate system (set of N = 400 virtual marker trajectories), while we determined local contour characteristics (e. ).Even under such extreme conditions, the algorithm produced stable virtual marker trajectories, i.e., trajectories without mapping violations.

Kymographs of local properties show characteristic patterns of amoeboid motility
The kymographs of local dispersion (LD), local motion, and curvature can be used to visualize, analyze, and quantify the expanding activity of a cell track.Note that the three quantities are closely related.For the persistently motile cell, the LD kymograph shows strong (positive) activity in a band-like structure along roughly half of the cell contour (width π), while the activity of local dispersion is less localized for the weakly motile cell and much less pronounced for the stationary cell.A similar scenario is seen in the local motion kymographs of the three cells.In broad terms, the local motion kymographs show more activity, e.g., areas of red and dark red color, than the LD kymographs.This is also apparent from a correlation plot of the two quantities (see S12 Fig).
In the following section, we chose the local dispersion as a basis to identify expansions being a product of local velocity and curvature.Another reason to choose LD as the quantity to define expansions is that many patches of high activity in the local motion kymograph contain multiple LD areas.The LD allowed us to divide these patches of high local motion into single separated expansions with high LD rate.

Virtual marker dispersion allows to identify and characterize expansions
Using the persistently motile cell track in In this context, we determine thresholds for expanding activity by dividing the 90th percentile of all positive LD values of the kymograph into three intervals of equal length.For contracting activities (LD < 0) the opposite thresholds were taken.By leaving out the largest and smallest values of the LD kymograph, the classification becomes less dependent on outliers.Additionally, we performed a prior smoothing of the kymograph as mentioned in the first paragraph of the previous section to reduce noise and, therefore, to reduce the number of small and separated patterns.
To highlight the events of the highest local expanding activity, we included positive local maxima of the LD kymograph (black dots) in panel (B).Local maxima falling inside regions with high or middle expanding activity were depicted as bold dots.Using the time-evolving coordinate system obtained from strongly regularized flow, September 9, 2021 18/31 the expansion/contraction areas shown in panel (B) are mapped back into the 2d plane of amoeboid motion, see panels (C) and (D), respectively.As a result, we obtain an automated visualization of the expanding activity during amoeboid locomotion.
The kymograph in panel (B) clearly shows the trace of expanding activity at the leading edge, located initially at around 3π/2 and then shifting towards π.At the same time, contracting activity occurs mainly at a distance of π from the leading edge.In panels (C) and (D), one can nicely see the explorative dynamics of the expansions at the cell front and the stably retracting back of the cell, the so-called uropod, where dark blue areas indicate faster contractions.Analogous graphics for a collection of motile and stationary cells can be found in S13 and S14 Figs, each containing 12 cell tracks.
In Fig 8, we present a close-up of two sequences of cell contours.The core expanding areas that shape the evolving cell contour can be seen clearly, e.g., in panel (A).The events of the highest expanding activity (black dots) seem to drive the September 9, 2021 19/31 expansion in many cases.Panel (B) illustrates strong contracting activity at the uropod, and the contraction of the membrane between two nearby expansions (blue area sandwiched between red areas).The opposite effect can be seen in panel (A) at the back of the cell, where a concave region between two convex contractions is identified as expansion.These patterns nicely illustrate the concept of local dispersion, being a product of the local velocity of virtual markers and the curvature along the contour segment.
In S15 Fig, every identified pattern with minimal growth time ∆t ≥ 3 is shown for the underlying cell track.

Statistical analysis of motility patterns
In this section, we illustrate the ability to statistically analyze a cell track based on our regularized contour flow approach.We used the persistently motile cell track for illustration.
Fig 9A shows the local dispersion kymograph of the persistently motile cell divided into four different phases (note that here the y axis starts at π/4).Until t = 200s, the cell moves upward with well observable expansions roughly between 1.1π to 0.2π, while the uropod is slightly above π/2.Then, the cell begins a phase of reorientation that lasts until t ≈ 280s, where larger expansions occur also at the former back of the cell.Subsequently, the cell changes its direction downward toward 3/4π to π (in the video it moves rightwards).In the last phase, the cell moves to the right-hand side by creating expansions at the front left, front right, and again front left of the cell.See S1 Video for a better understanding of the cell track.
We analyzed the expanding and contracting areas shown in Fig 7C and 7D with respect to the activity level (medium/high), duration, and position along the cell contour.In addition, we investigated the differences between expansions and September 9, 2021 20/31 The identified expanding and contracting areas are naturally partitioning by the sequence of contours into smaller 'slices' (see, e.g., Fig 8).We defined the area growth of an expansion/contraction as the area of the slice divided by the frame rate δt.We observed that the overall change of cell area attributed to expansions of high activity is substantially larger than for contractions of high activity.This illustrates that the cell motility of this cell track is driven by a higher number of fast (and potentially explorative) expansions, and a small number of fast contractions.Since in broad terms, the total area gain balances the total area loss, it further illustrates that area loss is to a larger extent attributed to slower and steadier contractions than it is for expansions.After the second phase (t > 280), when the cell has completed reorientation, the area change for both, expansions and contractions increased and even more so at the beginning of the last phase (t > 370).
Finally, in panel (C), the number of simultaneous expansions and contractions is shown (high activity only).In the first phase of the upward moving cell, the number of expansions is much higher than the number of contractions (mean: 1.52 vs. 1.12).In addition, we determined the fraction of time f T exp>2 and f T contr>2 with more than two simultaneous expansions and contractions, respectively.In the first phase, it is f T exp>2 = 0.146 vs. f T contr>2 = 0.035.In the following phases, these fractions increase for contractions: f T contr>2 = 0.198; 0.195; 0.341 (for phase 2; 3; 4), whereas for expansions the fraction is slightly smaller in the second phase and larger in the last two phases: f T exp>2 = 0.123; 0.207; 0.217 (for phase 2; 3; 4).This illustrates that the increased activity, as seen in panel (B), goes along with an increased number of expansions and contractions.Moreover, this indicates a more explorative character of the cell motion in phases three and four, while the cell seemed to be more stabilized  In other words, the local dispersion distribution shows only values that are larger than the thresholds for high expanding activity or smaller than the threshold for high contracting activity (the thresholds are ±0.087).The distribution of the LD of expansions further extends towards large values than the corresponding distribution of contraction towards small values (median of 0.13/s vs. −0.12/s).
In panels (B) and (C), the distributions of local motion and the curvature are shown for all virtual markers within expanding and contracting areas of high activity (same areas as above).For the local motion, we observed major peaks around 0.29 µm/s and −0.21 µm/s for expansions and contractions, respectively.Minor peaks correspond to inward expansions and outward contractions discussed in the previous section and shown in Fig 8 (see red expanding areas within a blue contracting region in panel (A), and a blue contracting area within an expanding region in panel (B)).In line with this observation, minor peaks of concave (negative) curvature and major peaks for convex (positive) curvature are shown.
Fig 10D shows the distribution of high activity expansions in direction of the moving cell.We identified two peaks in direction of the cell movement (front-left and front-right).Another peak is located at the back of the cell.A similar behaviour was presented for pseudopods in [30,31], where two different types of pseudopods were distinguished: (left/right) splitting pseudopods and de-novo pseudopods.By comparing the correlation between the growth time and the area of expansions and contractions in panel (E), we observed that expansions possess an area often twice as large as contractions with similar growth times.This indicates the difference between the faster and more explorative character of expansions at the cell front and the slower and stably retracting character of the uropod.This is also in line with the corresponding activities shown in Fig 9B .Finally, in panel (F) the distribution of growth times is shown for both, high activity expansions and contractions.Note that we used the term 'growth time' for both, expansions as well as contractions, owing to the fact that also contracting areas can be moving in an outward direction, as discussed earlier.The majority of growth times fall inside a range of 0s to 10s.Nevertheless, there are patterns, especially long persistent uropods, with growth times much larger than the range presented in these histograms.The average growth time of pseudopods observed in [30] is much higher (12.8 s) than the growth times of expansions presented in this work (4.9 s).This is not surprising, since our definition of expansions also takes short-lived objects into account.For example, the average number of expansions per minute for the persistently motile cell track (≈ 17.0) is much higher than the average frequency of pseudopods per minute (2.9 ± 0.2) as observed in [30].Additionally, concave regions between two contractions as in Fig 8A were detected as expanding areas as well, which raises the total number of expansions.See S13 and S14 Figs for similar graphics of a collection of 24 cell tracks, also containing the other two tracks from Fig 6.

Application to other kinds of cell motility
In this section, we demonstrate the flexibility of our method to study other kinds of cell motility.To this end, we extracted sequences of cell contours from videos of early embryonic killifish cells (Fundulus heteroclitus) [32][33][34] and keratocytes cultured from Central American cichlids (Hypsophrys nicaraguensis) [35, S1-S3 Movies].While the images of the embryonic killifish cells were obtained from fluorescence microscopy, bright-field microscopy was used to record the keratocytes.As before, the images were segmented by using a modified version of the active contour (snake) algorithm described in [21,26].
In Fig 11, three cell tracks for both applications, embryonic killifish cells and keratocytes, and corresponding local dispersion kymographs are shown.In panels (A) and (C), protrusion-driven cell motility can be observed for the embryonic killifish cells.Moreover, rotating waves around the leading edge (so-called circular movements) can be seen in panel (B) and in the last third of panel (C) as diagonal lines.In the chosen parametrization, right upward diagonals indicate counterclockwise rotations of protrusions while right downward diagonals indicate clockwise rotations.Notably, the switching between these kinds of locomotion can be nicely seen in panel (C) at t ≈ 630.In addition, we refer to the video accessible on [33] the cell track from panel (B) with dominant circular movements.
On the contrary, keratocytes possess a steady and persistent type of cell migration with only minor shape deformations.However, keratocytes can also show a more crawling kind of cell motility depending on the adhesion strength.In [35, S1-S3 Movies], cell track recordings for different adhesion strengths are provided: intermediate adhesion strength (panel (D)), low adhesion strength (panel (E)), and high adhesion strength (panel (F)).For low and intermediate adhesion strength, persistent cell track can be observed with a slightly higher local dispersion at the cell front for the case of low adhesion strength.For high adhesion strength, the cell migrates slower and exhibits larger cell shape deformations.Again, diagonal lines in panel (F) indicate circular movements as described earlier.Due to segmentation/mapping difficulties resulting from jumps in the video recording, pronounced artifacts can be seen as vertical lines at t = 360s and t = 750s.
These test applications underline the flexibility of our algorithm, handling contours obtained from fluorescence images as well as bright-field recordings.Moreover, the September 9, 2021 23/31

Discussion
With the ever-increasing amount of live cell imaging data and the continuously growing computational power, computer-automated techniques to analyze the morphology of cells have steadily developed over the past two decades.In particular, in cell motility research, morphological characteristics are commonly used to pinpoint phenotypic differences between mutant cell lines thus highlighting the mechanistic role of individual components of the underlying signaling pathways.While many static measures of cell shape have already been introduced early on [36,37], dynamic measures that quantify the temporal evolution of the cell shape proved to be more difficult to implement.First attempts focused on temporal changes of the projected cell area to deduce overall protrusion rates, for an example see [38].These approaches were later refined by local measures of cell boundary motion along selected line segments perpendicular to the cell border [39].Also, local space-time plots have been defined in this way [40,41].However, in all these cases, the direction of interest or the local placement of the kymograph had to be chosen manually, which severely limits a reliable long-time tracking of more complex cell shapes and introduces an arbitrariness related to the manual processing.The most promising approach to overcome this limitation relies on an active contour (snake), a closed chain of connected nodes (virtual markers) that is placed along the cell perimeter [42].Different rules have been proposed to propagate the markers from one contour to the next.In some cases, the distance between virtual markers is kept constant and markers are added or removed as the contour evolves [43].Other approaches that are inspired by mechanical spring models or concepts from electrostatics keep the number of markers constant and allow for local variations in the marker spacing [20].Here, we can distinguish two limiting cases.On the one hand, equidistance is enforced and markers on adjacent contours are connected in a one-to-one mapping by minimizing the sum of square distances between pairs of connected markers [26].On the other hand, markers are propagated from one contour to the next in a normal direction (normal flow) while the marker spacing evolves without constraints.The latter approach has been implemented using a level set method to cope with problems related to finite time sampling [19].However, high computational costs and the rapid buildup of highly uneven marker distributions limit the use of the level set method in practical applications.
In this article, we have introduced a family of marker flows that incorporates these different cases into a general framework.In particular, the regularization that we introduced in Eq (16) includes the two extreme scenarios described above as limiting cases.In the limit of large λ, we obtain an equidistant mapping, whereas, in the limit of small λ, we approach the shortest path flow, respectively, the reverse normal flow.Tuning λ allows us to systematically shift between these two limits.
Once a flow of virtual markers is computed on the evolving cell contour by any of these methods, it defines a coordinate system, in which different local quantities can be displayed, such as curvature, membrane displacement, or the intensity of a fluorescently tagged membrane-associated protein.Essential to our approach is the clear separation of local quantities derived by weakly regularized flows between two consecutive contours and the coordinate system which is based on a single strongly regularized (global) flow onto which each quantity of interest is mapped.This way, trajectories of each quantity can be obtained for the entire time period.Note that for all of these coordinate choices it is generally acknowledged that the dynamics of the virtual membrane markers do not reflect the motion of specific membrane lipids or proteins, as the membrane itself is a very complex and dynamic structure [19,26].In particular, lateral flows may occur due to membrane recycling, so that the dynamics of individual molecules or domains in the membrane do not necessarily correspond to September 9, 2021 25/31 morphological changes and will prevent a one-to-one mapping of molecular markers on adjacent contours.
Amoeboid motion is primarily driven by localized membrane protrusions, so-called pseudopodia.Identifying and tracking pseudopodia has thus been an important focus of the morphodynamic analysis of amoeboid cells.The first substantial pseudopod statistics were generated by computer-assisted manual image processing, relying on the expert judgment of the investigator [30,44].From this, a first automated software routine for pseudopod tracking was developed [22] and successfully applied also to analyze the chemotactic navigation of amoeboid cells [31].It relies on a complex decision tree that defines pseudopodia based on a sequence of threshold criteria applied to the local curvature, the virtual marker movement, and the local area change.In this way, the frequency and direction of pseudopod formation, their sizes, lifetimes, and other quantitative measures were extracted.While successfully providing a first quantitative database of pseudopod characteristics, this approach has the drawback that it requires the choice of several parameters that are tuned to the characteristic properties of pseudopods in starvation-developed D. discoideum cells.If cells display protrusions with a more diverse range of shapes and time scales, a reliable tracking is difficult to achieve with this approach.
Later, a more compact criterion for the detection of localized protrusions was proposed [21,26].It relies on thresholding a distance measure between the current position of virtual markers and the cell boundary at a later time point.The calculation of this distance measure, however, lacks a clear underlying definition and is computed in an ad hoc fashion for the specific data set (see supplementary material of Ref. [21]): First, each virtual marker is mapped from its current position onto the closest point on the future cell contour.As this generates a highly non-uniform distribution of target markers, with protrusive areas particularly poorly covered, the target markers on the new contour are then redistributed by two successive smoothing steps, using averaging windows of specific sizes.The time interval between the successive contours for the distance projection, as well as the smoothing parameters for redistribution of the target markers, were hand-picked by the investigator.Note, however, that this could be envisioned as one step in an iterative procedure to minimize our cost functional.
In the present work, we introduce a novel approach to define, identify and analyze localized expansions on dynamically evolving contours of amoeboid cells.expanding areas are defined via a single threshold value.In contrast to previous approaches, we chose the virtual marker dispersion as the underlying quantity, since it combines information on both, the marker displacement and the local curvature.The marker dispersion is mathematically well defined by Eq (25) and does not require additional empirical smoothing steps.Based on this criterion, we not only detect individual expansion events but we capture the entire shape and the complete temporal evolution of an expansion in a fully automated fashion, see An implementation of the methodology, obtaining kymographs from regularized flows in general and detecting expansion/contraction patterns from these kymographs in particular, are fully accessible and well documented in our software Package AmoePy.We also provide the data of multiple cell tracks and simple artificial test cases on which the algorithm was validated.Finally, all figures and results from this article can be easily reproduced in AmoePy.
As outlined in the Introduction, the overall aim is a quantitative, data-driven model of amoeboid motility.The presented theoretical framework is a first step in this direction.Because of its rigorous mathematical formulation, its efficient avoidance of mapping violations during larger shape deformations, and its moderate computational costs, it is a suitable choice for such a model.We envision that point events of high expanding activity may be used to define a point process in the space-time coordinate September 9, 2021 26/31 system.To reflect the often observed persistence in motility, so-called self-excited Poisson cluster processes or Hawkes processes may be favorable choices.The point process can serve as a skeleton for expansion activity that is 'completed' to a random realization of a kymograph, based on the statistics of a local quantity.We illustrated this idea in Fig 12, where we used the local motion statistics (B) to reconstruct a local motion kymograph (D).The idea is to use such realizations of kymographs to reconstruct a cell track and eventually to assimilate the time-lapse microscopy data into a mathematical model of amoeboid motility.

Fig 1 .
Fig 1. Virtual marker flows for two test cases w.r.t different degrees of regularization: strongly regularized (A, D), weakly regularized (B, E), and non-regularized (C, E).In panel (A)-(C), a stretching effect is illustrated as gray area.Furthermore, mapping violations are highlighted as red lines (F).
Multiple testing routines based on artificial test data and experimental data, see S1 Text and S4 Fig for more details, the computation time drastically; see S1 Text and S5 Fig for more details.Fig 2 shows the pseudocode to compute a regularized flow for a given regularization parameter λ.Being able to determine the regularized flow for a given parameter λ is the prerequisite to finally compute the global (strongly regularized) and local (weakly regularized) flows.These are based on the two regularization parameters λ glo ≫ λ loc > 0. The local quantities are then represented in the coordinate system of the global flow, yielding the kymograph representation.A pseudocode describing the computation of these kymographs is shown in Fig 3.

Fig 3 .
Fig 3. Algorithm to compute curvature, local motion and local dispersion.The two regularization parameters λ glo and λ loc are used in the RegFlow routine (see Fig 2) to determine the global (strongly regularized) and local (weakly regularized) flow.The output comprises the smoothed contours Γ 0 , . . ., Γ K−1 as well as geometric quantities of interest such as curvature, local motion, and local dispersion.

Fig 4 .
Fig 4. Cell trajectory of a persistently moving amoebae.(A) Fluorescence image with closed string of M = 400 equidistant nodes resulting from the segmentation process; shown is only every fifth segmentation point (blue points).(B) Smooth representation Φ k of the cell contour (orange line) obtained by spatial Gaussian regression on the segmentation points.Every fifth cell contour is displayed as dashed gray line.(C) Entire cell track of K = 500 cell contours (only every fifth shown).(D) Global trace of the cell track (gray area) and the trajectory of the center of mass of the contour (solid line, color coded as in panel (C)).The initial contour is shown as dashed black line and the final contour as dashed gray line.

Fig 5 .
Fig 5. Impact of regularization on the distribution of virtual markers for (A) no regularization (λ = 0), (B) weak regularisation with λ = 0.1, and (C) strong regularisation with λ = 1000, illustrated on two frames (roughly 20s apart for illustration purpose).Using the strongly regularised so-called coordinate markers as a means to map local characteristics into a kymograph, the lower panel shows the local motion (D) and curvature(E).The local motion is defined by the magnitude of each mapping vector, which are determined based on the weakly regularised marker flow.All panels correspond to the persistently motile cell of Fig 4.
g. local motion or local dispersion) based on a re-initialized weakly regularized flow.The local characteristics were subsequently represented in the coordinate system obtained from a strongly regularized flow.In panel (D), the local motion obtained from the re-initialized weakly regularized flow with λ = 0.1 is shown for the same cell track as in Fig4and the entire time-lapse microscopy recording (500 frames, δt = 1s).The local motion clearly shows regions of fast-moving membrane parts at the leading edge (red areas) and at the back of the cell (blue areas).The curvature kymograph in panel (E) shows characteristic lines of strongly convex (orange) and concave (green) membrane parts.Inclined lines of curvature (and local motion) may result from adherent parts of the cell moving along the cell contour as well as shifting effects of virtual markers due to arc length changes.It is important to notice that a kymograph depends on the underlying time-evolving coordinate system.In contrast to the procedure mentioned above, one may also compute local characteristics along the global flow without re-initializing from a uniform distribution of markers after every time step.However, this is not recommended, because either local information gets lost (strongly regularized flow) or clustering and thinning effects of markers become too prominent (weakly regularized flow); see S9 Fig for global flows based on different regularization parameters and the resulting kymographs.As shown in S10 Fig, we further challenged our algorithm by computing a strongly regularized (global) flow of the cell track from Fig 4 based on a few frames only (8 out of 500 Fig 6 shows cell tracks and corresponding kymographs for three different motility patterns: (A) persistently motile cell; (B) weakly motile cell; and (C) an almost stationary cell.All kymographs are smoothed by a Gaussian filter with a standard deviation of three markers in space and a standard deviation of one contour in time.See S11 Fig for the same kymographs but without smoothing.Moreover, videos of all three cell tracks and their corresponding kymographs can be found in S1-S3 Videos.

Fig 6 .
Fig 6.Comparison of different cell tracks of Dictyostelium discoideum: persistently motile (left), weakly motile (middle) and almost stationary (right).The corresponding kymographs contain information on the local dispersion (left), local motion (middle) and curvature (right).For details see text.

Fig 4 ,
we describe next, how to use the LD to define expansion areas and expanding events.Based on the local dispersion kymograph in Fig 7A, we defined areas of medium (light red) and high (dark red) expanding activity as well as medium and high contracting activity (light and dark blue, respectively), see Panel (B).

Fig 7 .
Fig 7. From the local dispersion kymograph to expanding areas and expansion events.(A) Local dispersion of a persistently motile cell as in Figs 4 and 6 (left-hand side).(B) Discretised local dispersion with different areas of activity: high (dark red), medium (light red), low (white) expanding activity, and low (white), medium (light blue), high (dark blue) contracting activity.Local maxima of positive local dispersion are depicted as black dots.Areas of medium and high expanding (C) and contracting (D) activity mapped back on the trace of the cell track.

Fig 8 .
Fig 8. Expanding and contracting areas with corresponding expansion events.Illustrative sequence of the contour dynamics for 96s ≤ t ≤ 144s (left), and 316s ≤ t ≤ 350s (right) based on the cell track shown in Fig 7. Features with high and medium expanding activities are shown in dark and light red, respectively.Features with high and medium contracting activities are shown in dark and light blue, respectively.All patterns shown possess a minimal growth time of 3s.The black dots show local maxima of the local dispersion in areas of medium and high expanding activity.

Fig 9 .
Fig 9. Statistical analysis of example cell track.(A) Local dispersion kymograph with thresholds as in Fig 7. (B) Area growth of cell segments which are part of identified expansions (red) and contractions (blue) of high and middle intensity.(C) Number of expanding and contracting areas with high intensity with respect to time.The time with > 2 features are colored gray to highlight the change in activity in the different phases.

Fig 10 .
Fig 10.Statistical analysis of example cell track.Distributions of local dispersion (A), local motion (B) and curvature (C) inside expanding and contracting patterns with high intensity.(D) Circular histogram displaying the angle, where high expanding activity appears along the cell contour.(E) Correlation between area and growth time of identified patterns.(F) Histograms of growth times of expansions and contractions with high intensity.

Fig 10
Fig 10  gives further insight into the expanding activity.Panel (A) shows the distribution of LD values of all virtual markers within expanding and contracting areas with high activity.In other words, the local dispersion distribution shows only values that are larger than the thresholds for high expanding activity or smaller than the threshold for high contracting activity (the thresholds are ±0.087).The distribution of the LD of expansions further extends towards large values than the corresponding distribution of contraction towards small values (median of 0.13/s vs. −0.12/s).In panels (B) and (C), the distributions of local motion and the curvature are shown for all virtual markers within expanding and contracting areas of high activity (same areas as above).For the local motion, we observed major peaks around 0.29 µm/s and −0.21 µm/s for expansions and contractions, respectively.Minor peaks correspond to inward expansions and outward contractions discussed in the previous section and shown in Fig8(see red expanding areas within a blue contracting region in panel (A), and a blue contracting area within an expanding region in panel (B)).In line with this observation, minor peaks of concave (negative) curvature and major peaks for convex (positive) curvature are shown.Fig10Dshows the distribution of high activity expansions in direction of the moving cell.We identified two peaks in direction of the cell movement (front-left and front-right).Another peak is located at the back of the cell.A similar behaviour was presented for pseudopods in[30,31], where two different types of pseudopods were distinguished: (left/right) splitting pseudopods and de-novo pseudopods.By comparing the correlation between the growth time and the area of expansions and contractions in panel (E), we observed that expansions possess an area often twice as large as contractions with similar growth times.This indicates the difference between Fig 8 and S15 Fig for example.

Fig 12 .
Fig 12. Future outlook.Underlying distributions are based on averaging over several cell tracks.(A) Circular histogram displaying the position, where expansions and contractions events appear along the cell contour.(B) Distribution of local motion of expansion and contractions events above a threshold ±1/3.(C) Simulated data obtained from self-excited Poisson point processes (so-called Hawkes processes) on the unit circle.Afterward, we obtained regions of high (red) and low (blue) intensity w.r.t. to a clustering algorithm.(D) Continuous kymograph obtained from a regression model (e.g.Gaussian process regression) based on sampled magnitudes at event locations shown in (C) from the distribution in (B).