Biophysical motivation & description of our contour dynamics approach

After obtaining smooth contours at discrete time steps the question of how to link them in time remains open. Mapping a marker on one contour to its nearest point on the consecutive contour is one desirable feature; achieving a regularly spaced distribution of markers on the consecutive contour is another desirable feature. Typically, a regularization parameter is used to balance the two features. In Fig 1, mappings defined by differently regularized flows of virtual markers are shown for two consecutive contours: a strong regularization (left column), a weak regularization (middle column), and no regularization at all (right column). While the strongly regularized flow preserves the distances of neighboring virtual markers, the non-regularized flow is defined by the shortest paths from the first contour to the second.

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).

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 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 R2, 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 S1 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 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

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 Φ.

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 ∂_θϕ_k(θ) < 0, is illustrated in Fig 1F, 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 θ₀ 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 θ₀, and let χ_k be the mapping of Γ₀ to Γ_k recursively defined by

This gives χ_k(θ₀) = θ_k for k = 0, …, K − 1. The Lagrangian description Ξ(tk,θ0)∈R2 is linked to the Eulerian description via

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.

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 ϕ_k via

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 μ₀(⋅) ≡ 1/(2π), this defines the transported density on all contours. In the Lagrangian picture this transport preserves the density μ₀(⋅), 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

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 ∫02πμ(θ)=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 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 H_k[ϕ_k] = F_k[ϕ_k] + λU_k[ϕ_k] 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.

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

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

yielding

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). If ϕ_k(θ_k,i+1) − ϕ_k(θ_k,i) = θ_k+1,i+1 − θ_k+1,i ≤ 0 for some k, we say that ϕ_k exhibits mapping violations.

In general, the distribution of virtual markers θ_k ∼ μ_k on Γ_k depends on the initial distribution of virtual markers θ₀ ∼ μ₀ on Γ₀, since the density μ_k results from the transport of μ₀ 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 Γ₀. 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

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 Γ₀.

We summarize the proposed numerical workflow:

Given the segmented contours γ₀, …, γ_K−1 as in Eq (1), determine the continuous representations Φ_k defined in Eq (2) satisfying the conditions in Eq (3).

Determine contour based quantities like the curvature in Eq (4) or the center of mass in Eq (5).

To determine the (strongly-regularized) coordinate flow, consider N equally spaced markers

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.

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 ξ₀, …, ξ_N−1 on Γ_k.

Determine contour flow based quantities like the local dispersion in Eq (25) or the local motion:

see also Eq (9).

Map all quantities, based on contour flow as well as contour only, onto the (strongly regularized) global flow.

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,

Multiple testing routines based on artificial test data and experimental data, see S1 Text and S4 Fig for more details,

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

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 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.

The algorithm input consists of three parameters r, σ_n, and N regarding the Gaussian process regression, a regularization parameter λ and the initial (segmented) contours γ₀, …, γ_K−1. The regularized flow is described by the output variables Γ₀, …, Γ_K−1 denoting smooth contours evaluated at a finite number of coordinate markers θ0eval,…,θK−1eval.

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 Γ₀, …, Γ_K−1 as well as geometric quantities of interest such as curvature, local motion, and local dispersion.

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 be measured.

Nevertheless, we provide the computation times of our algorithm “RegFlow” shown in Fig 2 on a standard computer, see S6 Fig. 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.