Three-component contour dynamics model to simulate and analyze amoeboid cell motility in two dimensions

Amoeboid cell motility is relevant in a wide variety of biomedical processes such as wound healing, cancer metastasis, and embryonic morphogenesis. It is characterized by pronounced changes of the cell shape associated with expansions and retractions of the cell membrane, which result in a crawling kind of locomotion. Despite existing computational models of amoeboid motion, the inference of expansion and retraction components of individual cells, the corresponding classification of cells, and the a priori specification of the parameter regime to achieve a specific motility behavior remain challenging open problems. We propose a novel model of the spatio-temporal evolution of two-dimensional cell contours comprising three biophysiologically motivated components: a stochastic term accounting for membrane protrusions and two deterministic terms accounting for membrane retractions by regularizing the shape and area of the contour. Mathematically, these correspond to the intensity of a self-exciting Poisson point process, the area-preserving curve-shortening flow, and an area adjustment flow. The model is used to generate contour data for a variety of qualitatively different, e.g., polarized and non-polarized, cell tracks that visually resemble experimental data very closely. In application to experimental cell tracks, we inferred the protrusion component and examined its correlation to common biomarkers: the F-actin density close to the membrane and its local motion. Due to the low model complexity, parameter estimation is fast, straightforward, and offers a simple way to classify contour dynamics based on two locomotion types: the amoeboid and a so-called fan-shaped type. For both types, we use cell tracks segmented from fluorescence imaging data of the model organism Dictyostelium discoideum. An implementation of the model is provided within the open-source software package AmoePy, a Python-based toolbox for analyzing and simulating amoeboid cell motility.


Introduction
Amoeboid movement belongs to the most widespread kinds of eukaryotic cell motility [1,2].It plays a key role in many biophysical and physiological processes such as wound healing, cancer metastasis, and immune system responses and is characterized by dynamic changes of the cell shape [3][4][5].In this context, the cell is moving forward by creating protrusions, so-called pseudopodia.The location of these protrusions and the frequency of their formation mainly define the overall trajectory of the cell [6].In addition to these explorative and faster membrane protrusions, the cell retracts its rear (so-called uropod) in a slower and steadier way.By this coordinated interplay of protrusions and retractions, the cell can move persistently and efficiently in a crawling-like fashion.Finally, the cell track is affected by the creation and rupture of adhesion contacts to the substrate [7].
On an intracellular level, amoeboid migration is initialized by changes of the actin cytoskeleton, a microfilament meshwork consisting of many multi-functional proteins.More precisely, the cytoskeleton grows by polymerization and shrinks by depolymerization of actin filaments.This rearrangement of the actin network is controlled by additional cytoskeletal proteins initializing capping, severing, and nucleation of actin filaments creating different meshwork formations such as bundling, cross-linking, and branching [8].The underlying mechanics are regulated by biochemical signaling pathways initializing different actions such as proliferation by cell growth and cell division, cell aggregation, and coordinated cell death [9,10].Furthermore, signaling pathways are linked to membrane receptors, enabling the cell to move persistently and directly towards a chemical source (chemotaxis) by sensing gradients and extracellular cues [11].
The formation of pseudopodia is the key aspect of amoeboid motility.The size and shape of these pseudopodia can vary significantly for different cell types [12].Some organisms such as Dictyostelium discoideum even possess different modes of locomotion and can switch back and forth between them [13].Therefore, pseudopodia are subdivided based on their nature into different types: lobopodia, lamellipodia, filopodia, reticulopodia, axopodia, invadopodia, and others [14][15][16].Pseudopodia can be also classified depending on the exact location of their appearance: Y-shaped (split) pseudopodia at the front, and de novo pseudopodia at the side or rear of the cell changing the direction of its movement [17].In this work, we focus on lobopodia, which are characterized by cylindrical and finger-shaped protrusions, and lamellipodia possessing a flat and broader structure [12].
A wide variety of amoeboid motility models have been proposed tackling different aspects of amoeboid migration such as membrane protrusion and retractions, trajectories of the cell's centroid, its polarization, and the influence of chemoattractant cues.Many proposed models are based on the concentration of interior biochemical compounds of the cell [18,19].In these reaction-diffusion models, motility patterns often result directly from the interplay of different processes controlling local excitation and global inhibition of the intracellular signaling and cytoskeletal dynamics [20][21][22][23].Reaction-diffusion models have been used to describe the self-organized polarization of the cell in the presence of chemoattractants [24][25][26][27][28], and the occurrence of intracellular waves and oscillations [29][30][31].A large number of these approaches are based on phase-field models which are used to describe the transition between different phases such as liquid and solid states or the interior and exterior of the cell.In the latter case, one modeling approach is to define the interior and exterior of the cell as binary states with a smooth transition function.The cell membrane is then defined where the transition function reaches the exact midpoint of both states [32].Phase-field models are often used to describe D. discoideum [13,26,[32][33][34], but also for other cell types [35][36][37][38][39].
Other approaches are mechanical models in which different forces affect the cell motility from within and from the outside.Mechanical models differ in complexity and dimensionality to target different problems, e.g., the formation of fibroblasts (1D) [40], the influence of contraction and adhesion sites to the substratum on amoeboid cell motility (2D) [41], and the evolution of the cell surface obtained from triangulation under chemotaxis (3D) [42].Furthermore, different physical methods and assumptions are used for the underlying model equations such as active gel physics where the gel consists of polymer filaments permeated by a solvent [43], hydrodynamics to model the internal cytoplasmic fluid [44], and the modeling of the extra-cellular matrix [45,46].Most importantly, mechanical models can be easily adjusted for unusual types of cell migration such as amoeboid swimmers [47,48].In [49], a mechanochemical model is proposed for which the underlying parameters are calibrated by using Bayesian optimization.
Finally, level set methods are used to simulate cell tracks, e.g., as part of a mechanical model of the cell cortex combined with an excitable network acting as activator/inhibitor system [19].The excitable network is triggered by random fluctuations and can be enhanced with additional gradient stimuli and polarization modules [50,51].In [52], the stochastic extension of pseudopods during chemotaxis is modeled.Furthermore, stochastic differential equations have been used to model the trajectory of the cell's centroid, e.g., by a (generalized) Langevin equation [53][54][55][56][57].In [58], the centroid trajectory was modeled by a Markov chain approach on a discrete domain obtained from hexagonal tiling.Geometric equations of evolving curves are commonly used to describe cell migration [23,42,59,60].A novel contour evolution method based on the curve-shortening flow (CSF) or, alternatively, with an optional additive function which is then called forced CSF, is presented in [61] and then applied to cell migration in [62].A comparison of the different model approaches mentioned above can be found in [63][64][65].
In contrast to the above approaches, our model is designed to infer key characteristics of individual cell tracks, i.e., the intensity of protrusions and retractions during amoeboid cell motility.The model is therefore intended to be simple and intuitively comprehensible to ensure a fast and straightforward estimation of underlying model parameters and to be capable of producing a specific motility behaviour for a given input.Especially the second property is sometimes hard to achieve with mechanical models, which tend to have a higher complexity with a large number of model parameters and entangled subprocesses.This makes it difficult to draw a direct link between model parameters and a desired motility outcome [46].By inferring the above characteristics, our model can also be used to identify differences between cells and to classify them.The second main goal of our model is to simulate a variety of quantitatively different and realistic contour dynamics producing cell tracks which can be hardly distinguished from experimental ones.Based on the model's October 25, 2022 3/32 capability to simulate such versatile contour dynamics, we deem it to be applicable to experimental cell tracks for varying degrees of motility and persistence, or even different types of locomotion.Briefly, our model evolves two-dimensional contours representing the cell membrane and is based on three components.The first component is driven by a stochastic term generating membrane protrusions and is modeled by (1) the intensity of a self-exciting Poisson point process (so-called Hawkes process [66][67][68]) or, alternatively, (2) an Ornstein-Uhlenbeck like diffusion process.In this context, we show that the Hawkes process, due to its self-exciting nature, is suitable to produce a cascade of protrusion events to ensure a persistent cell migration.The second two components are mathematically well-defined geometric flows initiating contour retractions: the area-preserving curve-shortening flow (APCSF) to regularize the shape/arc length of the contour; and a further flow introduced as area adjustment flow (AAF) which expands/shrinks the cell contour with respect to a specified reference area.For more information on the APCSF, see [69][70][71].Based on the interplay of the above components, the model defines the formation of protrusions as well as retractions.It is therefore linked to other mechanistic models evolving the cell contour as an elastic object in time and space.
In the following, we demonstrate how our model can be used to simulate realistic cell tracks.Then, we analyze experimental cell tracks (D. discoideum) by inferring the model-based protrusion and retraction components.In this context, we compare the inferred protrusion component with commonly used biomarkers: the density of filamentous actin close to the membrane and its local motion.An implementation of the model as well as a graphical user interface to simulate cell tracks is provided in our Python-based toolbox AmoePy [72].

General notations and underlying coordinate system
The notation and theoretical framework used in this work have been established in [73].Primarily, this framework is used to (1) obtain smooth contour representations from stacks of segmented microscopy images and (2) track reference points (so-called virtual markers) between successive contours.More precisely, smooth contours and the corresponding contour curvature can be derived easily from a discrete set of segmentation points (so-called active contour or snake) by using a Gaussian process regression (GPR), a Bayesian regression model based on kernel functions.In our case, we used a Poisson kernel function as underlying covariance function Briefly, the choice of r cont in k rcont (•, •) affects the rigidity and stiffness of the resulting contour.Furthermore, an additional noise parameter σ noise specifies the deviation between the initial segmentation points and the regression function obtained from the GPR.First, we consider K ∈ N cell contours denoted by Γ k with k = 0, . . ., K − 1.The corresponding time points are denoted by t k = k • δt with δt > 0. The coordinates of these contours are given by the following mapping where θ denotes the normalized arc length coordinate and Φ Noteworthy, connecting consecutive contours in time and space is intrinsically not well defined, i.e., there are multiple ways to do so [74][75][76].In many approaches, varying constraints are introduced to track reference points/virtual markers from one contour to the next one, e.g., by using electrostatic field equations [74], level-set methods [74,75], or mechanistic spring equations [76].Naive contour mapping approaches include the propagation of virtual markers (VM) based on shortest paths to the next contour or by choosing paths in normal direction only.However, these approaches can change the order of neighboring virtual markers (so-called topological mapping violations).
In [73], we address this issue by proposing a novel regularizing family of contour flows, connecting consecutive contours while preserving desirable characteristics of the underlying mapping trajectories.Depending on a regularization parameter λ reg ≥ 0, this family of contour flows includes the two extreme cases: either enforcing shortest VM trajectories between contours (λ reg = 0) or preserving equal distances between neighboring VM's for every time step (λ reg ≫ 0).For any flow between two contours Γ k and Γ k+1 , a mapping is induced, which describes the translation along the contour under the flow.By assuming we ensure that φ k is a one-to-one mapping, i.e., mapping violations between Γ k and Γ k+1 do not occur.Now, we can define a virtual marker trajectory starting at θ 0 ∈ [0, 2π) by the iteration In the following, we used a Lagrangian reference frame to describe geometric flows and the resulting evolution of virtual markers on the contour.This Lagrangian reference frame is denoted by χ k and recursively defined by: For the limit of infinitely dense contours, we introduce the following notations Given a VM p 0 = Φ 0 ( θ) on the first contour Γ 0 with arc length coordinate θ ∈ [0, 2π), we can now track this VM in time and space which corresponds to the function t → Φ(t, θ).Finally, we introduce a virtual marker distance ratio defined as: For the discrete case of N ∈ N virtual markers θ k,0 , . . ., θ k,N −1 on the contour Γ k , this ratio can be rewritten as: measuring the distance between neighboring virtual markers divided by the distance of equidistantly spaced VM's.
October 25, 2022 5/32 Contour propagation vs. contour mapping.Our model contains two different dynamics: (1) the contour propagation based on a model function f : R + × S 1 → R and (2) the contour mapping under which the driving force of our model is transported.Briefly, we evolve contours for short time periods by propagating an equidistant set of contour points with respect to for an arc length coordinate θ ∈ [0, 2π) and time t ∈ [0, T ].A propagation of these markers in tangential direction affects the position of the markers on the contour but does not affect the shape of the contour.Hence, the contour propagation in Eq (6) is observed in normal direction only.However, such a normal flow would lead to thinning and clustering effects of markers over time.For this reason, we use a second kind of dynamics, namely regularizing flows described as in [73] to ensure an evenly-spaced distribution of virtual markers.This flow is used to transport the underlying driving force in our model.Furthermore, it constitutes the underlying coordinate system necessary to draw graphical representations (so-called kymographs) of each model component.In contrast to the contour propagation, the virtual markers under the regularizing flow are propagated also in tangential direction.In S1 Text, we present a detailed description regarding the implementation of the model.In S1 Fig, we illustrate the two different kinds of marker trajectories: (1) the contour propagation (green dashed lines) and ( 2) the contour mapping (blue dashed lines) under which the stochastic protrusion process X prot (t, θ) is transported and the underlying coordinate system is based on.

Regularizing the shape and size of contours via geometric flows
Our model to evolve two-dimensional cell contours is based on three components: a protrusion term based on a stochastic process accounting for membrane protrusions and two geometric flows accounting for membrane retractions by regularizing the shape and area of the contour.Due to the separation into one protrusion component and two retraction components, the model provides an independent handling of these two key features of amoeboid cell motility.The first geometric flow is defined by the area-preserving curve-shortening flow (APCSF) and denoted by f APCSF .Briefly, the APCSF evolves a two-dimensional contour to a circle while preserving the area enclosed by the initial contour.It therefore minimizes the arc length of the contour without affecting the contour area.In the absence of other influences, the evolution of every contour to a circle is energetically favorable to reduce the surface tension of the membrane.For example, this behaviour can be observed during cell death or by treating cells with Latrunculin to dissolve the actin cytoskeleton [77].The second geometric flow, which we call area adjustment flow (AAF), is denoted by f AAF .The AAF shrinks/expands the contour towards a predefined reference area.Since the underlying contour data relies on two-dimensional cross sections of three-dimensional cells, the resulting contour area can change significantly.Therefore, a certain change of the contour area is desirable and possible in our model.By using both flows, we achieve a regularizing effect on the arc length (APCSF) and the area (AAF) necessary to counteract the forward movement initialized by the protrusion component f prot .This component is based on a stochastic process, e.g., a Hawkes intensity process or an Ornstein-Uhlenbeck process, and should be strictly positive since it describes the formation of protrusions only.
In Fig combination of all three components.In the first row, four different samples drawn from the same stochastic process X prot are displayed.The first two cases (panels (B) and C)) lead to significantly growing contours since the area adjustment is absent.In comparison, the contour dynamics in panel (C) are much smoother due to the regularizing effect of the APCSF.In panel (D), dynamics with highly curved but similarly sized contours are produced.However, self-intersections can easily occur without the APCSF.Therefore, a combination of all three components, as in our model, is necessary (last row).
Area-preserving curve-shortening flow.To regularize the contour arc length in our model, we used the APCSF: where κ(t, θ) denotes the curvature and n(t, θ) the outward-pointing normal vector at time t ∈ R + for a virtual marker with arc length coordinate θ ∈ [0, 2π) on the first contour.Here, L(Φ(t, •)) denotes the total arc length and is defined by the following functional: The APCSF is defined as the gradient flow of this functional under the area-preserving constraint A(t) := A(Φ(t, •)) = A(0) for all t > 0, where the contour area is defined by: October 25, 2022 7/32 As a consequence, the APCSF evolves every contour to a circle of the same area, minimizing the contour arc length to L(t) t→∞ −−−→ 2 πA(0) while maintaining its area.
Area adjustment flow.While the APCSF has no effect on the contour area, the protrusion component would expand the cell contour most of the time and therefore also its area.The area adjustment flow counteracts this expansion.It is defined as with reference area A ref ∈ R + and center of mass trajectory Φ CM : R + → R 2 defined by For a choice of the reference area A ref , we take the 1st percentile of the entire area time series of an experimental cell tracks.Other choices are possible.
It is easy to see that an area adjustment is achieved by this flow.Furthermore, the AAF affects the contour in normal direction only and is shape-preserving if used without the other two components.In S1 Text, we introduce another area regularizing flow which is based on the gradient flow of the area functional.This flow minimizes the contour area much more rapidly by affecting the contour in normal direction only.However, this flow is not shape-preserving.

Hawkes intensity process as protrusion component
Statistical analyses of different sequences of cell contours have shown that a protrusion event increases the probability of nearby follow-up protrusions [17].We therefore modeled the protrusion component f prot in our model by the intensity of a self-exciting Poisson point process, a so-called Hawkes process.Due to the self-exciting nature of the Hawkes process, a cascade of protrusion events can be generated which result in substantial and persistent contour changes.In this way, we obtain contour dynamics with a significantly moving center of mass instead of a fluctuating membrane only.
We used a product kernel function g(t, θ) = g 1 (t) • g 2 (θ) with temporal component g 1 (t) and spatial component g 2 (θ).The temporal kernel function is given by with arrival intensity α > 0 and exponential decay rate β > 0. As spatial kernel, we used the von Mises distribution, with κ M > 0 as concentration parameter and I 0 (κ M ) denoting the modified Bessel function of order 0. For a realization of the Hawkes process, the protrusion process X prot : R + × S 1 → R is defined as: with any spatio-temporal kernel function g : [0, T ] × [0, 2π) → R + , time scaling factor c s > 0, and VMDR as in Eq (4) accounting for local contour/arc length changes.For the sake of simplicity, we choose c s = 1s and the same kernel function as above, i.e., g(t, θ) ≡ g(t, θ).Hence, we can rewrite the protrusion process X prot in terms of the Hawkes intensity process λ(t, θ) from Eq (9): to highlight that X prot directly resembles the intensity of the Hawkes process and not the underlying realization of point events or the Hawkes process itself.
Based on the underlying choice of parameters, different motility characteristics can be adjusted with our model: • the general movement speed by w prot , • the number of protrusions by λ 0 and α, • the duration of protrusions by β, • the size of protrusions (many small protrusions vs. a single large protrusion) by κ M , • membrane fluctuation vs. creation of pseudopods with a substantial movement of the center of mass by α, • non-polarized vs. polarized contour dynamics by r pol .
In S3 Fig, the kernel functions g 1 (t) and g 2 (θ) are illustrated for varying parameters α, β, and κ M as well as the Hawkes intensity λ(t, θ) for a fixed set of parameters.In S1 Text, we illustrate an alternative approach by modeling the protrusion component with an Ornstein-Uhlenbeck type of diffusion process.
October 25, 2022 9/32 Three-component contour dynamics model In this section, we formulate our contour dynamics model based on the three components: a protrusion component based on a stochastic process X prot from Eq (13), the APCSF from Eq (7), and the AAF from Eq (8).In our model, the following parameters are required: • weight parameters w prot , w APCSF , w AAF > 0, • a reference area A ref > 0, • additional parameters regarding the stochastic process X prot .Furthermore, the computation of the following geometric quantities is necessary: For a virtual marker Φ(0, θ) with initial arc length coordinate θ ∈ [0, 2π), its trajectory t → Φ(t, θ) with t ∈ [0, T ] is given by the differential equation where f : R + × S 1 → R is defined as with the following three components: In Fig 2, these three components are displayed.In panel (A), the protrusion component is displayed with a generated sample of the underlying stochastic process X prot in the top right square.X prot is defined on a unit circle and is mapped onto the contour with respect to the contour arc length and a reference point θ = 0 (blue dot).In panel (B), the APCSF is shown with its final state (dashed grey contour): a circle with the same area as the initial contour.Since the APCSF is mainly defined by the contour curvature, the contour shrinks faster for regions with large positive (convex) curvature and expands faster for regions with large negative (concave) curvature.In panel (C), it can be observed that the area adjustment is weaker for contour segments, which are closer to the center of mass.The final state is again displayed as a dashed gray contour.
Finally, in panels (D-F), we display test scenarios of simulated contour dynamics each based on a single component only.In panel (D), we observe a steadily growing contour since the f prot is the only active component.Under the APCSF (panel (E)) the contour evolves to a circle by expanding concave parts and shrinking convex parts of the contour.In the case of AAF being the only active component, we observe a shape-preserving shrinkage of the contour.In this case, the reference area (A ref = 60µm 2 ) was set to be smaller than the initial contour area.Alternatively, by choosing a larger reference area, a shape-preserving growth of the contour would occur.

Inference of model components and parameter estimation
Besides simulating contour dynamics, our model is used to infer the different model components f prot , f APCSF , f AAF of experimental cell tracks to analyze and characterize cell tracks on the individual level.Importantly, we want to quantify the intensity of protrusions and retractions separately from each other.Usually, this is done by computing the local motion which, however, comprises the characteristics of the entire contour dynamics: protrusions, retractions, as well as minor membrane fluctuations.Finally, by inferring the underlying model components and estimating the corresponding component weights, our aim is to classify cells based on different motility types.
Since the APCSF and AAF from Eq (16) are based on the contour curvature, arc length, and area; f APCSF and f AAF can be computed separately for each time step solely based on the current contour.For this reason, one can explicitly determine the only remaining component f prot to propagate one contour to the next one.This approach enables us to replicate the experimental cell track with our model for an estimated set of model parameters.This set includes the reference area A ref and three model component weights w prot , w APCSF , and w AAF .
While the APCSF does not affect the contour area, the (positive) protrusion component increases the contour area.Therefore, the reference area A ref should be chosen relatively small to achieve a counteracting shrinkage effect of the AAF.In our case, we haven chosen the 1st percentile of the entire area time series of the cell track.
October 25, 2022 11/32 Next, we estimate f APCSF and f AAF and their corresponding weights, w APCSF and w AAF , by making use of the local motion kymograph of the underlying cell track.More precisely, we determine w APCSF and w AAF such that negative regions (i.e.retractions) of the local motion are optimally captured by the above components.Initially, the remaining component f prot is estimated by deducting f APCSF and f AAF from the local motion kymograph.Afterwards, we choose w prot such that the sample variance of the underlying protrusion process X prot is standardized with Var(X prot ) = 1.In a next step, we further tune the inferred protrusion component f prot by minimizing the distances of virtual markers propagated with respect to f prot and the "receiving" contour at the next time step.In this optimization step, we use the built-in least-squares method from the Python package SciPy.Finally, we evaluate the goodness of fit of the inferred protrusion component.Since f APCSF and f AAF are computed in advance, the remaining protrusion component corrects any missing contour dynamics to propagate one contour to the next one.For this reason, the inferred protrusion component can be negative if a contour retraction is stronger than the APCSF and the AAF have predicted.Therefore, a good fit is given if the estimated protrusion component is (mostly) positive, i.e., the contour retractions are successfully captured by the other two components.

Hawkes process based simulations of amoeboid cell motility
In this section, we show that the Hawkes process is suitable to simulate ameboid cell motility.With the proposed model a variety of qualitatively different contour dynamics were generated.
The protrusion component in this section is based on a Hawkes process defined as in Eq (13).The underlying model weights were set to w prot = 7.5µm 2 /s, w APCSF = 0.1µm 2 /s, and w AAF = 1µm/s.As reference area we have chosen A ref = 80µm 2 .First, we simulated contour dynamics based on a non-polarized test scenario realized by choosing r pol = 0.Then, we have chosen a polarized test scenario by choosing r pol = 0.5.A summary of all parameter values is displayed in Table 1.As initial contour we have chosen a circle with an area equal to A ref .
Fig 3 shows exemplary non-polarized contour dynamics, which are driven by a Hawkes process.We observed that the Hawkes process can enforce a substantial change of the center of mass trajectory as displayed in panels (A) and (B).The self-excitation can be also observed in the kymograph of the protrusion component f prot (panel (D)).In this kymograph, point events realized from the Hawkes process are depicted as circles and show clusters as expected from the self-excitation property.The corresponding Hawkes intensity was then used to define the protrusion process in our contour dynamics model as in Eq (13).The same red patterns of the f prot kymograph can be also found in the final local motion kymograph (panel (C)).The retractions of the contour dynamics are mainly defined by the area adjustment f AAF (panel (E)).By comparing this kymograph with the area plot in panel (G), we see that dark blue patterns, indicating strong retractions, occur when the contour area A(t) is much larger than the reference area A ref .
Due to the definition of the APCSF component f APCSF , the curvature at each part of the contour can be inferred from the corresponding kymograph (panel (F)).Since the contour dynamics start from a perfect circle, possessing a constant and relatively small curvature, the kymograph starts with values close to zero.Later on, blue horizontal stripes indicate the position of protrusions and the rear of the cell, which possess a large positive curvature and are retracted therefore by the APCSF.In

AAF weight
List of parameters for artificial cell tracks based on a Hawkes process.
contrast, concave regions of the cell contour correspond to red horizontal stripes since they are enforced to expand under the APCSF.Finally, the influence of the APCSF, minimizing the contour arc length, is shown in the plot in panel (H).sufficient to reproduce the zigzag movement.For this case, the clustering of point events (panel (D)) is even more pronounced than for the non-polarized case.From the local motion kymograph, we infer that the protrusions (red regions) occur mainly at θ = π defining the front of the cell.On the other side, retractions (blue regions) occur near θ = 0 and θ = 2π.In contrast to the non-polarized scenario, the curvature lines in the kymograph representing f APCSF are less horizontal.Instead, the characteristic diagonal stripes stand for protrusions created at the front of the cell which are then moved along both sides of the contour until they reach the rear of the cell.This has also been observed experimentally in [78].
In addition to the Hawkes process, we also modeled the creation of protrusions by a simple Poisson point process.We generated cell tracks for each of the two scenarios described above where the protrusion component is only defined by the first generation of point events without further offspring as for the Hawkes process.To achieve the same number of point events, we increased the background intensity from λ 0 = 1s −1 to λ 0 = 5s −1 .
In S4 events with less significant event clusters.Therefore, for the non-polarized scenario, membrane fluctuations without large displacements of the center of mass were observed.For the polarized scenario, the cell moves persistently in one direction without having major turns, zigzag movements, or additional de novo pseudopodia.Hence, the Hawkes process was better suited to model amoeboid cell motility than a standard Poisson point process.
Animations of the contour dynamics and the corresponding kymographs from The effect of the self-excitation of the Hawkes process can be clearly seen in form of cascades of point events determining the direction and the movement of the cell.
We furthermore examined a second alternative to the Hawkes process, an Ornstein-Uhlenbeck type process, see S1 Text for more details.However, the resulting simulations are less satisfactory due to a higher number of protrusive features that are October 25, 2022 14/32 not observed in experiments and additional artifacts such as a pulsating membrane and a swimming type of locomotion.Since a self-exciting property is missing in this approach, it is much more difficult to generate cascades of multiple protrusions and subsequent reorientation phases of the contour dynamics.Contour dynamics and the corresponding model components obtained from this approach are displayed in S1 Text with animations shown in S3 Vid.
How to choose the regularization parameter λ reg ?To ensure stable contour dynamics within our model, we used regularized flows defined as in [73].Based on a regularization parameter λ reg ≥ 0, we can enforce a more even distribution of virtual markers for every time step/contour.This way, thinning/clustering effects of virtual markers over time can be avoided, see S1 Text and S1 Fig for more details.
In this section, we investigated the influence of this regularization on the overall contour dynamics and studied under which circumstances clustering/thinning effects of virtual markers as well as other artifacts can occur.Moreover, we challenged our model by varying the temporal resolution.In this case, the underlying contour mapping is also affected since the number of contours is increased/reduced.
In Fig 5, we present simulations of non-polarized cell tracks (panel (B)) as well as polarized cell tracks (panel (C)) generated with the parameter choices in Table 1 but for varying regularization parameter λ reg ∈ {0.01, 0.1, 10, 1000}.The trace of each cell track (gray area) and the corresponding centroid trajectories (colored lines) are displayed.For the non-polarized case (panel (B)), the overall contour dynamics were substantially altered in cases of weak (λ reg = 0.01) or strong regularization (λ reg = 1000).For a medium regularization (λ reg = 0.1, 10), the model produced contour dynamics similar to each other with almost identical contours at the end of each track.Analogous observations were made for the polarized test cases in panel (C).By decreasing λ reg , the overall motility is also reduced.This can be understood from Eq (13): a low regularization results in thinning effects of VM's at the leading edge which increases the virtual marker distance ratio (VMDR) and therefore decreases the protrusion process.In panel (C), we noticed differences in the main direction of the cell track due to different contour mappings at the beginning of each track.However, the overall contour dynamics was not affected by increasing λ reg .
In panels (D) and (E), histograms are displayed showing the effect of each regularization scheme on the distribution of virtual markers for non-polarized (left) and polarized (right) contour dynamics.Each histogram displays the relative frequencies with respect to the virtual marker distance ratio from Eq (5) for all contours of a simulated cell track.In case of a weak regularization, the thinning and clustering effect are reflected by a higher variance of the VMDR with distances up to 2 times larger than in the equidistant case.Especially in panel (E) for λ reg = 0.01 (light blue histogram), a major peak at VMDR ≈ 0.5 indicates a strong clustering of virtual markers at the rear of the cell.In case of a strong regularization, an almost even distribution of virtual markers was achieved for every time step/contour which is reflected by VM distance ratios close to 1.
In panels (F) and (G), kymographs of the APCSF component corresponding to the non-polarized (left) and polarized (right) tracks are displayed.For a weak regularization (λ reg = 0.01), prominent thinning and clustering effects of virtual markers were observed, see dashed box in panel (G).In case of a strong regularization (λ reg = 1000), an equidistant distribution of virtual markers is enforced at each time step.As a consequence, an emerging protrusion has a direct effect on all virtual markers along the cell contour leading to a cyclic shift of all markers.These shifts were observed as rotating elements of the cell contour indicated by sharp diagonals in the APCSF kymograph, see dashed boxes in panel (F).For medium regularization choices (λ reg = 0.1, 10), the above artifacts (clustering effects, rotating elements) were not observed.
A summary of our findings is presented in Fig 5(A), where the influence of λ reg on the distribution of virtual markers is shown.In this context, we computed the standard deviation σ VMDR of the virtual marker distance ratio for varying λ reg .For a weak regularization λ reg = 0.01, a loose connection of neighboring virtual markers is observed which leads to thinning and clustering effects.If the regularization is chosen too strong (λ reg = 10 2 , 10 3 ), other artifacts such as rotating elements of the contour might occur.For this reason, we recommend a medium regularization: 0.1 ≤ λ reg ≤ 10.In our simulations, we have chosen the upper limit λ reg = 10 to obtain a more equidistant distribution of virtual markers which ensures stable contour dynamics even for a longer time period T > 500s while avoiding rotating artifacts of the contour.
For more details, see S6   As mentioned above, clustering and thinning effects of virtual markers were prominent for λ reg = 0.01.For the other choices of λ reg , we observed that all kymographs are barely affected at all, resulting in similar contour dynamics (top right).Major differences in the cell's direction are noticed only at the beginning of the cell track when a stable distribution of virtual markers is not yet reached.
In S8 Fig and S9 Fig, we present simulated non-polarized and polarized cell tracks for varying temporal resolution δt = 0.25, 0.5, 1, 2, 2.5, 3. 3s and fixed regularization parameter λ reg = 10.We observed that the local motion kymograph and, hence, the general contour dynamics are not substantially affected.In some cases, the overall direction of the contour dynamics was altered due to small differences at the beginning of the track.
In summary, we observed that by varying the regularization parameter within the range λ reg ∈ [0.1, 10], the overall contour dynamics are barely affected.The same observation was made for a varying temporal resolution.Based on the above studies, we haven chosen a medium regularization λ reg = 10 and a relatively dense temporal resolution of δt = 0.5s for the following simulations.

Long-time simulations are stable and show a normal diffusive behaviour.
To demonstrate the capability of our model to produce stable contour dynamics even for a longer time period, we simulated a variety of cell tracks under both scenarios, the non-polarized and the polarized case as described above.that qualitatively and quantitatively different cell tracks can be generated with our model, we provide a rationale for the inference approach in the following section, where we applied the model on experimental cell tracks and different types of locomotion.
In Fig 6, we present the diffusive behaviour of 50 cell tracks, respectively, for both polarization scenarios.For each scenario, a selection of ten cell tracks are highlighted in different colors (see panels (A) and (F)).The center of mass trajectories of these cell tracks for T = 500s are displayed in panels (B) and (G).Furthermore, an enlarged excerpt of panel (B) is shown in panel (C).The root mean squared displacements (RMSD) of the trajectories are depicted as gray circles at time steps t = 0, 100, . . ., 500s.The corresponding mean squared displacement (MSD) is shown in panels (D) and (H), indicating normal diffusion for the non-polarized case and a ballistic regime (so-called superdiffusion) for the polarized case for the first 500 seconds.In panels (E) and (I), the MSD is displayed for a longer time period of T = 10000s, clearly showing a linear relation between time and MSD indicating normal diffusive behavior also for the polarized case, see also S10 Fig, where the center of mass trajectories of the polarized case are presented for the time period of T = 10000s.Again, the RMSD is indicated as gray circles (δt = 2000s).Finally, based on linear fits (red dashed lines), we computed the diffusion coefficients: D ≈ 0.16µm 2 /s for the non-polarized case and D ≈ 34.33µm 2 /s for the polarized case.
In S4 Vid and S5 Vid, the contour dynamics of the respective cell tracks from

Inferring the protrusion component from experimental data
In the first part, we used our model to simulate a variety of different cell tracks based on stable and, in particular, realistic contour dynamics.In the second part, we will now apply our model to experimental data to analyze cell tracks on the individual level and to classify them based on different types of locomotion: the amoeboid type and a so-called fan-shaped type.The underlying microscopy imaging data was published in [79].As mentioned earlier, the APCSF and AAF from Eq (16) only depend on the current contour and deterministic quantities such as contour curvature, arc length, and area.Hence, the remaining component f prot that propagates one contour to the consecutive one, can be determined explicitly.For a sequence of experimental cell contours, this approach enables us to infer the different model components for a given set of model weights; see S1 Text for more details on how to estimate these parameters.
In Fig 7A, contour dynamics of D. discoideum are displayed for a time period of T = 500s and a frame rate of δt = 1s.This cell track is based on fluorescence microscopy data (see panel (B)), from which the cell contours are segmented (colored lines).As mentioned, we have chosen the 1st percentile of the entire area time series as underlying reference area, i.e., A ref = 91.23µm 2 .By following the approach from S1 Text, we estimated the following model weights: w prot = 6.634 µm 2 /s, w APCSF = 0.057 µm 2 /s, and w AAF = 3.532 µm/s.In comparison to parameter values estimated for other cell tracks, see S11 Fig and S12 Fig, we observed that w prot and w AAF are relatively large while w APCSF took a medium value.This observation coincides with the following cell track characteristics: a fast and persistent movement, an area time series within a smaller range 90µm 2 < A < 120µm 2 (see panel (F)), and contour dynamics showing an intermediately regularized curvature.The accuracy of our computational approach is demonstrated in panel (B), with virtual markers (black circles) being effectively propagated onto consecutive contours.In the following, we compare the inferred protrusion component (panel (E)) to commonly used biomarkers: the local motion of the membrane (panel (C)) and the F-actin density/fluorescence intensity close to the membrane (panel (D)).
Inferred protrusion component in comparison to other biomarkers.In contrast to the commonly used local motion which reflects the entire contour dynamics, i.e, protrusions, retractions and minor membrane fluctuations, the proposed model separates the creation of protrusions, which depend on f prot , and retractions, which depend on f APCSF and f AAF .We expect the protrusion component to be correlated to the F-actin density near the membrane.However, here we have chosen images, where the fluorescence signal saturates at high F-actin levels to facilitate segmentation of the cell contours.Thus, by following this approach, the resulting fluorescence intensity of the F-actin density provides less details.
In Fig 7, kymographs of all three quantities are displayed: the local motion (panel (C)), the fluorescene intensity (panel (D)), and the inferred protrusion component (panel (E)).We expect that a kymograph of a successfully inferred protrusion component would be strictly positive and only indicate regions where protrusions occurred.The remaining retractions are then covered by the other two components of our model: the AAF and the APCSF.In the local motion kymograph, local protrusions (red regions) and retractions (blue regions) are nicely shown, while the fluorescence kymograph captures the main characteristics of the cell motility only.The protrusion component kymograph mainly displays protrusive areas (red regions) with only a very few negative regions (blue).We thus conclude that the missing retractive regions are successfully captured by the AAF and the APCSF.
In S13 Fig, we present our approach to measure the fluorescence intensity near the membrane by averaging over ellipses along the cell contour.To improve the quality of the fluorescence intensity kymograph, we post-processed the recorded microscopy images by using a (tenfold) image upsampling, i.e., we smoothed all microscopy October 25, 2022 19/32 images by increasing the number of pixels tenfold.However, the resulting kymograph showed no differences to the kymograph based on the original image data.
In S14 Fig, the protrusion component from Fig 7 is displayed as well as the remaining components f APCSF and f AAF .As the underlying set of model weights, we used the same estimates mentioned above.The contributions of all three model components to the overall dynamics are shown for each time step/contour.In this context, we observed that the APCSF accounts for approximately 10% of the overall dynamics.In contrast, the AAF and the protrusion component are more prominent accounting for approximately 40% and 50%, respectively.This may indicate that the APCSF mainly resolves a smaller number of strongly curved contour segments, while slower contour retractions at the rear of the cell are controlled by the AAF.
In addition, we present kymographs of the same quantities for an alternative set of model parameters, where positive values of f prot are more favored; see S1 Text for more details.For this case, f prot and f AAF accounted for the overall dynamics equally, whereas f APCSF was observed to be negligible most of the time.
October 25, 2022 20/32 Impact of a-priori parameter choices on the analysis of experimental data.
The F-actin density near the membrane is often interpreted as a marker of protrusive activity.Thus, we would expect some correlation between the F-actin density and the protrusion component in our model.Since the latter is influenced by the choice of the model weights of the two other model components AAF and APCSF, the question arises to what extent the correlation depends on the a priori choices of w AAF and w APCSF .We reformulated our model to be based on relative weights r prot , r APCSF , and r AAF = 1 − r prot − r APCSF and an overall velocity parameter w f instead of the absolute weights w prot , w APCSF , w AAF used in Eq (16); see S1 Text for more details.In S15 Fig, we display protrusion component kymographs for varying parameters r prot ∈ {0.05, 0.5, 0.8}, r APCSF ∈ {0.01, 0.05, 0.1}, and w f ∈ {1, 5, 10, 20} of the cell track from Fig 7 and the corresponding correlation coefficient between these kymographs and the fluorescence intensity/F-actin density.For a relatively strong APCSF, i.e. r APCSF > 0.1, we observed distinct positive and negative horizontal patches induced by the contour curvature.On the other hand, we obtained predominantly positive values in the kymograph for a strong AAF, e.g., the top left kymograph for w f = 20µm/s.Since the AAF affects/shrinks every part of the contour, the counteracting protrusion component is increased to the same amount for the entire contour in order to replicate the given contour dynamics.On the other hand, for a too small AAF, we observed distinct negative (blue) regions in the protrusion component kymograph, indicating that additional retractive forces are necessary to replicate the cell track.
By comparing these protrusion component kymographs with the local motion and the fluorescence intensity from Fig 7 panels (C) and (D), we observed the following similarities and differences.In this context, we focus on the two examples highlighted as black and white dashed boxes in each kymograph.The first example (300s < t < 400s) clearly indicates a retraction as shown in the local motion kymograph.However, a significant remaining density of F-actin is seen in the fluorescence intensity kymograph.This pattern is also displayed for most of the protrusion component kymographs in S15 Fig, e.g., the top left kymograph on the second page (w f = 5µm/s, r prot = 0.05, r APCSF = 0.01).This means that our model predicts a significant amount of f prot necessary to slow down the retraction such that the modeled contour dynamics resemble the experimental data.Without this contribution of the protrusion component, the APCSF and AAF would enforce an even stronger retraction.A similar relation can be observed in the second example (t > 400s), but to a lesser extent as in the first example.For all cases, the resulting correlation coefficients between protrusion component and F-actin density fell within a range of [0.32, 0.49].The best fit with a correlation of ρ = 0.49 was obtained for two parameter choices of {w f , r prot , r APCSF } with distinctly different values, e.g., {10, 0.5, 0.01}, and {20, 0.8, 0.01}.In case of the estimated (relative) model weights used in Fig 7 given by {10.2, 0.65, 0.01}, we achieved a comparable correlation coefficient of ρ = 0.47.Since the correlation coefficient is invariant under changes in location and scale, we conclude that most of the protrusion kymographs displayed in S15 Fig differ in magnitude primarily.However, if the APCSF is chosen too strong (r APCSF ≥ 0.1), the resulting protrusion kymographs show substantial differences and a decreasing correlation coefficient.For this reason, the APCSF weight should be chosen relatively low r APCSF ≤ 0.05.
In general, we observed that the protrusion component inferred by the model is correlated to the underlying F-actin density.While the protrusion component is affected significantly by the parameter choice, the impact on the above correlation coefficient is minor.
Classification of contour dynamics based on cell motility types.In the first part of this work, we have shown that our model can be used to simulate a variety of different contour dynamics.Here, we demonstrate that our model can also be applied to different experimental cell tracks.By inferring the protrusion component and estimating the underlying model weights, we analyzed contour dynamics on the individual level and classified them based on two different types of locomotion: the amoeboid type and the so-called fan-shaped type.Again, the contour dynamics were derived fluorescence images of D. discoideum, where frames are recorded with a temporal resolution of δt = 1s (amoeboid type) and δt = 4s (fan-shaped type).
The contour parameters r cont = 0.6, σ noise = 0.05 were chosen as in the simulation part, for more information see S1 Text of [73] In contrast to the local motion (first kymograph row), the protrusion component (third kymograph row) contains only a few negative (blue) regions, which means that the APCSF and AAF capture most of the contour retractions successfully.Furthermore, we see strong correlations between all three quantities.However, the protrusive regions in the fluorescence intensity kymographs are larger than for the other two kymographs, which is a direct consequence of detector saturation during fluorescence imaging.From the estimated model weights, we can infer that the third cell track is more motile (w prot = 6.505µm 2 /s) than the other two cell tracks (w prot = 4.513µm 2 /s) for both), which coincides with the contour dynamics displayed on top.Furthermore, we see that the AAF weight w AAF directly influences the variability of the area time series, i.e., less variability for the second cell track (w AAF = 1.429µm/s) vs. more variability for the first and third cell track (w AAF = 0.902µm/s and 0.676µm/s).Finally, we can observe that the elongated contour shape of the second and third cell track coincide with higher APCSF weights (w APCSF = 0.048µm 2 /s and 0.070µm 2 /s) compared to the smaller weight of the first cell track (w APCSF = 0.023µm 2 /s) In the case of the fan-shaped cells in S12 Fig, we see clear blue stripes in the local motion kymograph at the rear of the cell, which corresponds to a smaller F-actin density in the fluorescence intensity kymograph.Because of the characteristic kidney-shaped cell contour, our model predicts a negative protrusion component at the cell's rear also displayed as a light blue stripe in the kymographs.This indicates that the APCSF and the AAF were not able to fully capture this retraction.Since the APCSF evolves every contour to a circle, it counteracts the characteristic concave-shaped contour of the cell.For this reason, we expect the estimates of w APCSF to be lower than for amoeboid locmotion.Indeed, the APCSF weight was estimated to be zero for all three cells.Therefore, by inferring the impact of the APCSF, we have shown that our model is capable of distinguishing fan-shaped cells from standard amoeboid cells.
For the second and third cell track similar weights w prot = 4.485µm 2 /s vs. 4.428µm 2 /s and w AAF = 1.844µm/s vs. 1.76µm/s were estimated, which is in accordance with the similar contour dynamics reflected by comparable mean centroid velocities of 0.091µm/s and 0.099µm/s, respectively.In contrast, we estimated a lower protrusion component weight w prot = 3.121µm 2 /s for the first cell track, which is in agreement with a smaller mean centroid velocity of 0.069µm/s.Moreover, we estimated a smaller AAF weight w AAF = 0.911µm/s, which might be a result of the ever-increasing contour area of this track (130µm 2 to 217µm 2 ).

Discussion
We developed a novel model to simulate and analyze the contour dynamics of amoeboid cell migration.The contour dynamics model is based on three components: (1) a stochastic protrusion component based on a self-exciting Poisson point process also know as Hawkes process, (2) an area-preserving curve-shortening flow (APCSF) to regularize the contour arc length, (3) and a further geometric flow introduced as Area Adjustment Flow (AAF) to control the contour area.While the first component controls the forward movement of the cell, the latter two components control contour retractions, occurring most often at the rear of the cell.
First, we have shown that our model is capable of generating a variety of cell tracks with different spatio-temporal patterns.We simulated non-polarized as well as polarized contour dynamics that are hardly distinguishable from experimental data and stable over a long time period.Secondly, we applied our model to experimental cell tracks to analyze key motility characteristics based on the inference of the three model components and the estimation of the underlying model component weights.By examining the correlation between the inferred protrusion component and the fluorescence intensity reflecting the F-actin density, we demonstrated that the creation of pseudopods is correctly accounted for by the protrusion component of our model.Furthermore, our model was capable of decoupling the explorative protrusions from the slower contour retractions, which in our model are solely based on the APCSF and the AAF.Finally, by estimating the model weights of each component, we demonstrated a simple approach to classify cells based on two locomotion types: the amoeboid and a so-called fan-shaped type.
The simplest approaches to model amoeboid cell motility focus on the motion of the center of mass of the cell [56][57][58].Often, the center of mass trajectory is modeled by a Fokker-Planck equation [27] or a Langevin equation [53][54][55][56][57]. Due to a lower model complexity, even a large number of cell tracks can be generated relatively easy and fast.These models are then used to examine statistical quantities such as the diffusion coefficient, the persistence time, and the drift velocity.While our model focuses on the contour dynamics, the center of mass trajectory and the corresponding statistics can be derived easily.For example, we determined a diffusion rate of D = 0.16µm 2 /s in our test case of simulated non-polarized cell tracks which is in accordance with experimental measurements [27,80].
Other approaches include biochemical models focusing on intracellular signaling molecules and cytoskeletal components [24,25] or specific biophysical mechanisms such as cell polarization during chemotaxis [18,19].Some models combine intracellular processes with specific membrane patterns such as the formation of pseudopods [22] or changes of the cell shape in general [81].Furthermore, phase-field models are often used in which the transition of different states such as solid/liquid or interior/exterior of the cell are described [26,33,34].In these models, biochemical markers such as the F-actin density are often propagated in time and space and drive displacements of the cell membrane [13,32].
In contrast to most of the existing approaches which focus on the microscopic level of chemical and biophysical process during cell migration, our model describes amoeboid motility on a macroscopic level.By relying on the deformation of an elastic object, namely the cell contour, it has a mechanistic basis.Usually, mechanistic approaches are based on the interplay of multiple forces affecting the cell from the outside but also from within [40,41,46].However, some mechanistic models rely on multiple entangled subprocesses making it more difficult to achieve a direct causality between a chosen parameter regime and a specific motility behaviour [46].Here, our model differs substantially due to its dependency on three components only and a low October 25, 2022 23/32 number of parameters: 4 model parameters and 4-5 additional parameters regarding the stochastic process.Desired motility characteristics such as the number, duration, and size of protrusions, the general movement speed or the cell polarization can be easily achieved based on an appropriate specification of the underlying parameter regime.Motivated by the biological insight that the location of protrusions depends on previous protrusions [17] and is triggered by a cascade of physiological events affecting the growth of the actin filament network [9], we have chosen a Hawkes process as the underlying stochastic process.Due to its self-exciting property, the Hawkes process is capable of generating multiple explorative protrusions in temporal and spatial proximity with subsequent reorientation phases characterized by a temporary decline of the cell motility.Compared to other processes such as the Ornstein-Uhlenbeck process or a standard Poisson process, the Hawkes process is therefore advantageous to model membrane protrusions.
Phase-field models often contains additional constraints to preserve the area/volume of the cell within a specific range [32,33,82].Constraints regarding the surface tension of the membrane are added to these models to regularize the shape or the perimeter of the membrane.In a similar fashion, we use the AAF and APCSF to regularize the contour area and curvature.In contrast to phase-field approaches, the APCSF and AAF affect and propagate a one-dimensional object only, namely the cell contour, in contrast to a two dimensional concentration of biomarkers or even the entire cytoskeleton.Furthermore, phase-field models are successfully used to study more complex processes such as cell division or the interaction and movement of multiple cells [83,84].
By focusing on the contour dynamics to model cell motility, further assumptions regarding the mapping of virtual markers along time and space were necessary.We assumed that the transport of the underlying protrusion process is based on the concept of regularizing contour flows which was previously introduced in [73].Other approaches to determine contour/virtual marker mappings include electrostatic field equations [74], level-set methods [74,75], or mechanistic spring equations [76].All of the above approaches, including ours, share the same problem that the mapping of consecutive contours is not defined a priori.We have shown that our contour mapping assumption has a minor effect on the overall contour dynamics and we provided a reasonable range of the underlying regularization parameter λ reg for which simulated contour dynamics are stable.
So far, only few approaches offer a full integration of the numerical model and experimental measurements, e.g., for fibroblast migration [49].Our model provides a full and automatized integration of experimental measurements due to its capability of inferring protrusion and retraction components individually and by offering a simple and straightforward estimation of the underlying model weights.The most common approach to tune computational models with respect to experimental data is to perform a sensitivity analysis [46,[85][86][87].In this context, different parameter regimes are chosen to simulate different motility behaviors and to determine the importance of each parameter on the simulated outcome.With a similar approach, we have shown that the inference of the protrusion component by our model is stable for varying model weights w prot and w prot .However, if the APCSF is too prominent, the inferred protrusion component is negatively affected, indicated by a weaker correlation with the underlying F-actin density.Regarding the classification of different motility types, we examined contour dynamics based on a single locomotion type (either amoeboid or fan-shaped).However, recent studies report spontaneous switching between these two locomotion types [88].
In summary, the proposed model sets new standards in simulating stable and, in October 25, 2022 24/32 particular, realistic contour dynamics.Due to a fast and straightforward parameter estimation and a fully automatized approach to infer protrusion and retraction characteristics, the model can be used to analyze experimental cell tracks on the individual level and to classify them.
(x) k (•) and Φ (y) k (•) the x and y coordinates of the contour, respectively.

Fig 2 .
Fig 2. Comparison of the three model components: f prot , f APCSF , and, f AAF .(A) Protrusion component driven by a stochastic process (generated sample shown in top right corner).The reference point θ = 0 is highlighted as blue dot.(B, C) APCSF and AAF with final states displayed in corresponding top right corners.(D-F) Example contour propagated individually by only one component for a time span of 1s.Regarding the AAF, the underlying reference area is given by A ref = 60µm 2 .

Fig 4
shows an artificial cell track based on a polarized Hawkes process.In contrast to the cell track of Fig3, the contour dynamics a characterized by a higher motility and a stronger persistence.Interestingly, we observed the zigzag pattern well known from experimental amoeboid cell tracks[17].The self-exciting behavior of the Hawkes process initialized with a spatially unimodal background intensity (see S2 Fig)is

Fig 3 .
Fig 3. Artificial non-polarized cell track with T = 500s based on a Hawkes process.(A, B) Contour dynamics (left, colored contours), the center of mass trajectory (right, colored line), and the trace of the entire cell track(right, gray area).(C, D) In the second row, kymographs of the local motion f and its protrusion component f prot are displayed.Point events realized by the underlying Hawkes process are depicted as circles in the protrusion kymograph.(E, F) In the third row, kymographs of the other two components f AAF and f APCSF are shown.Finally, in panels (G, H), the evolution of the contour area as well as the contour arc length are presented.
Fig 3 and Fig 4 are shown in S1 Vid and S2 Vid, respectively.In these videos, each point event generated by the Hawkes process is illustrated as a circle in the protrusion component kymograph (second row) and the animated contour dynamics (bottom left).

Fig 4 .
Fig 4. Artificial polarized cell track with T = 500s based on a Hawkes process.The contour dynamics, the center of mass trajectory, the contour trace, the kymographs of f , f prot , f APCSF , and f AAF as well as the contour area and arc length are shown in the same order as in Fig 3.
Fig, where the above non-polarized cell tracks are shown for a wider selection of varying regularization parameters λ reg ∈ {10 −2 , 10 −1 , . . ., 10 3 } and with all corresponding model components.Again, for the intermediate cases λ reg ∈ [0.1, 10], we noticed only very few differences in all shown kymographs.However, because of small changes of the contour mappings at each time step, the overall direction of the cell track can vary over time.The polarized cell tracks under varying regularization schemes are shown in S7 Fig.

Fig 5 .
Fig 5. Influence of regularization parameter λ reg on simulated cell tracks.(A) The impact of λ reg on the virtual marker distribution described by standard deviation of the virtual marker distance ratio σ VMDR and its reciprocal.(B, C) Non-polarized (left) and polarized (right) cell tracks generated with varying regularization parameter λ reg ∈ {0.01, 0.1, 10, 1000}.The same protrusion component is underlying in (B) and (B), respectively.(D, E) Histograms of the VMDR for varying λ reg for the non-polarized (left) and polarized (right) cases.(F, G) Kymographs of the APCSF component of all non-polarized (left) and polarized (right) cell tracks with regions of interest shown as black and white dashed boxes.

FurthermoreFig 6 .
Fig 6.Diffusion analysis of artificial cell tracks generated from a non-polarized test scenario (top row) and a polarized test scenario (bottom row).(A, F) For each scenario, ten exemplary cell tracks are shown.(B, C, G) Center of mass trajectories (colored lines) with root mean squared displacement (grey circles) at different time points with ∆t = 100s.While panel (B) and (G) are scaled equally, panel (C) is an enlarged excerpt of panel (B).(D, H) Corresponding MSD of these trajectories (bold black line) for medium time spans of T = 500s.(E, I) MSD for longer time spans of T = 10000s with linear fit (dashed red line).
Fig 6A and F are shown.

Fig 7 .
Fig 7. Inferring protrusion component from D. discoideum cell track.(A) Persistently motile cell track for T = 500s.(B) Microscopy image with fluorescence intensity (white to green color scheme) and segmented cell contours (red lines, every tenth shown).(C) Local motion kymograph showing expansions (red areas) and retractions (blue areas).(D) Relative fluorescence intensity with regions of high and low F-actin density displayed in red and blue, respectively.(E) The underlying protrusion component inferred from our model for a given set of model parameters.The resulting propagation of virtual markers from one contour to the next one are depicted as black circles in panel (B).(F) The contour area with predefined reference area A ref = 91.23µm 2 (dashed gray line).Finally, regions of interest are displayed as black and white dashed boxes.
. The reference area A ref was chosen individually for each cell track based on the 1st percentile of the entire time series of the contour area (gray dashed line): 62.43, 122.99, and 64.61 µm 2 for the amoeboid cell tracks; 142.15, 203.45, and 189.24 µm 2 for the fan-shaped cells.Moreover, the model weights were estimated for each cell track individually, described as in S1 Text.In S11 Fig, three tracks of an amoeboid motion as in Fig 7 are presented.

S13Fig.
Computation of relative fluorescence intensity for experimental microscopy data and tenfold upsampled data via ellipses along the cell contour.(PDF) S14 Fig. Model components extracted from experimental cell track of Fig 7 and their proportion on the overall velocity of the contour dynamics for two pairs of model weights.The model weights were estimated by minimizing sums of squared residuals: S (left column) and S + (right column), see S1 Text for more details.(PDF) S15 Fig. Protrusion component f prot extracted from the experimental cell track of Fig 7 for varying relative weights r prot ∈ {0.05, 0.5, 0.8} (vertical axis) and r APCSF ∈ {0.01, 0.05, 0.1} (horizontal axis) as well as varying overall velocity parameter w f ∈ {1, 5, 10, 20} (page axis).(PDF) S1 Vid.Contour dynamics with corresponding kymographs of a non-polarized cell track based on a Hawkes process as shown in Fig 3.The point events generated by the Hawkes process are depicted as white circles.The cell track is displayed at fivefold speed.(MP4) S2 Vid.Contour dynamics with corresponding kymographs of a polarized cell track based on a Hawkes process as shown in Fig 4. The point events generated by the Hawkes process are depicted as white circles.The cell track is displayed at fivefold speed.(MP4) S3 Vid.Contour dynamics for different modifications of the underlying Ornstein-Uhlenbeck process as protrusion process.The cell tracks are displayed at tenfold speed.(MP4) S4 Vid.Contour dynamics of non-polarized cell tracks based on a Hawkes process as shown in Fig 6.The cell tracks are displayed at tenfold speed.(MP4) S5 Vid.Contour dynamics of polarized cell tracks based on a Hawkes process as shown in Fig 6.The cell tracks are displayed at tenfold speed.(MP4) 1, simulated contour dynamics are shown based on: the protrusion component f prot only, f prot combined with f APCSF , f prot combined with f AAF , and a

Table 1 .
Choice of parameters and meaning.