# Scaling Manipulation Learning with Visual Kinematic Chain Prediction

Xinyu Zhang Yuhan Liu Haonan Chang Abdeslam Boularias
{xz653, yl1834, hc856<EMAIL_ADDRESS>
Rutgers University

###### Abstract

Learning general-purpose models from diverse datasets has achieved great
success in machine learning. In robotics, however, existing methods in multi-
task learning are typically constrained to a single robot and workspace, while
recent work such as RT-X requires a non-trivial action normalization procedure
to manually bridge the gap between different action spaces in diverse
environments. In this paper, we propose the visual kinematics chain as a
precise and universal representation of quasi-static actions for robot
learning over diverse environments, which requires no manual adjustment since
the visual kinematic chains can be automatically obtained from the robot’s
model and camera parameters. We propose the Visual Kinematics Transformer
(VKT), a convolution-free architecture that supports an arbitrary number of
camera viewpoints, and that is trained with a single objective of forecasting
kinematic structures through optimal point-set matching. We demonstrate the
superior performance of VKT over BC transformers as a general agent on Calvin,
RLBench, Open-X, and real robot manipulation tasks. Video demonstrations can
be found at https://mlzxy.github.io/visual-kinetic-chain.

> Keywords: Multi-Task Robot Learning, Manipulation

##  1 Introduction

There are numerous techniques in machine learning and computer vision that can
successfully learn a single general-purpose model from multiple diverse
datasets [1]. In robotics, however, despite the recent advances in multi-task
learning that enable a single policy to perform various tasks through
imitation learning with language instructions [2, 3], these methods are
typically constrained to a single robot and workspace. Some recent work, such
as RT-X [4, 5], leverages a large vision-language model (VLM) to directly
train policies with Behavioral Cloning (BC) on the Open-X Embodiment [5], a
collection of datasets crowd-sourced from various environments and robots.
However, these techniques require a non-trivial action normalization procedure
to bridge the gap between the different action spaces in diverse setups, such
as end-effector poses, joint positions, and velocities in various world
frames. This manual engineering procedure is currently custom-designed for
each training dataset, which affects the generalization and interpretability
of these models.

Therefore, a key question is: can we find an action representation that is
precise and universal for various setups and robots? To this end, we propose
the visual kinematic chain, which is the projection of the robot’s high-
dimensional kinematic structure into the image plane as pixels. Instead of
predicting the low-level robot actions, we propose to learn and visually
forecast the movement of robots’ kinematic chains in the image plane. This
visual approach provides a unified action representation for different robots,
and requires no manual adjustment since the visual kinematic chains can be
automatically obtained from the robot’s model and camera parameters.

To forecast the kinematic structures of various robots, we render the
kinematic chains into point sets by sampling points along links and performing
optimal matching [6] to minimize the earth moving distance [7, 8] between the
predicted point sets and ground-truth kinematic structures. Further, we
propose the Visual Kinematics Transformer (VKT), a convolution-free
architecture that is solely based on the attention mechanism in the RGB space,
which supports an arbitrary number of camera viewpoints. Our VKT is trained
with a single objective of forecasting kinematic structures through optimal
point set matching without seeing any low-level robot actions. VKT is deployed
to a specific environment by simply training a tiny head while freezing the
backbone. VKT demonstrates superior performance over BC transformers as a
general agent on Calvin (success length 1.08 vs 0.33), RLBench (success rate
61.7% vs 24.3%), Open-X, and real robot manipulation tasks.

To summarize, our contributions are threefold. (1) We propose the visual
kinematics chain as a precise and universal representation of quasi-static
actions for learning from diverse robot configurations. (2) We propose VKT, a
convolution-free architecture that supports an arbitrary number of camera
viewpoints and is trained with a single objective of forecasting kinematic
structures through optimal point set matching. (3) We present a thorough
empirical study of the performance of VKT on specialized and general agents
over a diverse set of language-conditioned tasks and environments.

Figure 1: Overview of the proposed framework. We use the visual kinematic
chain as a universal action representation across diverse robots and setups.
We propose the Visual Kinematics Transformer (VKT), an architecture based
solely on attention layers, which predicts the future movements of the visual
kinematic chain in images from multiple viewpoints. Our VKT is trained with
the earth-moving distance as the single learning objective without knowing any
low-level robot states or actions. When deployed in a specific environment, we
freeze the VKT as a backbone and attach a tiny head to project the VKT output
to actual robot commands.

##  2 Related Work

Robot Learning from Diverse Environments. Learning a single universal multi-
task policy on various robots, cameras, and task configurations is a
challenging problem that is receiving increasing attention [9]. One recent
strategy to solving this problem consists in training a network that takes as
input a robot’s kinematic structure [10, 11]. For instance, MetaMorph [12]
tokenizes the robot’s kinematic tree, which is then used as a transformer
prompt to produce per-joint action. However, these methods are only evaluated
in simple walking environments of Mujoco [13] without any variations in the
task, the camera viewpoint, and the workspace. Several recent works have also
looked into leveraging a large vision-language model (VLM) to directly train
robots, through behavior cloning, on a collection of crowd-sourced datasets
such as Open-X Embodiment [14, 15, 5]. Despite their impressive results, these
methods require a non-trivial action normalization procedure to reduce the gap
between action spaces such as end-effector poses, joint positions, and
velocities in various world frames. However, this manual engineering procedure
is custom-designed for each training dataset, which severely affects the
generalization and interpretability of these models. In contrast, our method
directly learns and visually predicts the movement of robots’ kinematic
chains.

Transformers for Manipulation Learning. Transformers are widely used in
manipulation learning [4, 16]. However, existing architectures require a fixed
number of camera viewpoints, which is not guaranteed when learning is
performed across various environments. Both RT-1/2 and RT-X select one
canonical view from each Open-X dataset [5, 4, 17]. RVT [18] supports multiple
but predefined viewpoints. PerAct [3] trains transformers over point-clouds
and is not limited by camera viewpoints. But methods based on point clouds
still suffer from the lack of datasets with depth information, as well as
scale discrepancies [19, 20]. In contrast, our visual kinematics transformer
(VKT) is a convolution-free architecture that is solely based on attention
mechanisms in the RGB space, without the need for depth data. VKT also
supports an arbitrary number of camera viewpoints.

Intermediate Action Representations. Instead of only using low-level robot
commands, incorporating intermediate actions into control policies has been
shown to increase data efficiency [21]. RT-H is a hierarchical policy that
predicts language motions such as “move arm forward” [22]. SWIM uses
affordance maps to unify grasping actions of humans and robots [21]. General
Flow employs a labeling process to translate end-effector motions into point-
cloud movements [23]. IMOP [24] uses invariant regions to describe key
manipulation poses. RT-Trajectory predicts the future 2D trajectory of the
end-effector’s center in a given image [25]. In contrast, our method forecasts
the movements of the entire kinematic chain in multiple camera viewpoints.
This representation provides a visually grounded action abstraction, while
capturing detailed and small motions at the same time.

Figure 2: Rendering a kinematic chain as a set of pixels in an RGB image

##  3 Visual Kinematic Chain Forecasting

We consider the problem of learning a single universal multi-task policy from
diverse robots and workspaces through a unified action space. Our key insight
lies in the use of the visual kinematic chain, which is the projection of the
robot’s high-dimensional kinematic structure into a set of pixels in the image
2D plane. The visual kinematic chain can be analytically generated with camera
parameters and robot models without any manual engineering. We show that this
simple representation allows us to train a single visual language-conditioned
policy over multiple environments. The VLM policy is trained with a single
objective of forecasting the future visual kinematic chain movements without
knowing any low-level robot actions. Our proposed transformer architecture is
built entirely with attention layers and can predict consistent multi-view
chains from an arbitrary number of camera viewpoints. To deploy the multi-task
policy in a specific environment, one simply needs to freeze the transformer
and train a tiny head to project the actions into actual robot commands. An
overview of our framework is illustrated in Figure 1.

Fitting any Structure with Point-Set Matching. The kinematic structures of
various robots differ significantly. To forecast different structures, we
propose to render them to point sets by uniformly sampling points along links,
as illustrated in Figure 2. Let $P=\\{p_{i}\\}_{i=1}^{M}$ denote the point-set
rendered from the ground-truth configuration of a kinematic chain, and
$Q=\\{q_{i}\\}_{i=1}^{N}$ denote the point-set predicted by VKT, where
$p_{i}\in\mathbb{R}^{2},q_{i}\in\mathbb{R}^{2}$. We first compute the pairwise
euclidean distance matrix $C\in\mathbb{R}^{M\times N}$, where
$C_{ij}=\parallel p_{i}-q_{j}\parallel$. We then minimize the earth-moving
distance $\operatorname{EMD}(P,Q)=\sum_{i,j}\gamma^{*}_{ij}C_{ij}$, where
$\gamma^{*}$ is a stochastic assignment matrix that matches points in $P$ and
$Q$. We find $\gamma^{*}$ using the Sinkhorn-Knopp algorithm to solve the
following optimization problem online for each mini-batch [26, 8].

$\begin{array}[]{r}\vspace{0.5em}\gamma^{*}=\arg\min_{\gamma\in\mathbb{R}_{+}^{M\times
N}}\sum_{i,j}\gamma_{ij}C_{ij},\text{ s.t. }\gamma
1=1;\gamma^{T}1=1;\gamma\geq 0.\end{array}$ (1)

Figure 3: Overview of our proposed visual kinematics transformer (VKT). For
each camera input, we encode the language instruction and RGB image as text-
and-vision tokens with CLIP [27]. Then, we concatenate the text tokens and the
kinematics tokens as query tokens. The kinematics tokens are learned
parameters. Next, the query and vision tokens are interweaved with a sequence
of our proposed multi-view dual attention block. For each block, the query
tokens are first updated with self-attention (orange). Then, cross attention
is applied with the query tokens as queries and the vision tokens as keys and
values (green). A cross-attention layer updates queries with keys and values,
we use ➝ to denote queries and $\multimapdot$ to denote keys and values. Then,
both the query tokens and the vision tokens are updated through cross-
attentions with tokens from other camera viewpoints (blue). Next, the vision
tokens are updated by cross-attention with query tokens as keys and values
(green). Finally, the $T$ kinematics tokens are projected into $T$ point sets
through an MLP, representing the visual kinematic chain in the current and the
future $T-1$ steps. The predicted point-sets are optimized through point-set
matching with the ground-truth, as shown in Equation 1.

Visual Kinematics Transformer (VKT). Our VKT network accepts a language
instruction and an RGB image as input and predicts a sequence of $T$ point-
sets. The dimensions of the kinematics tokens are $T\times D$, where $D$
denotes the embedding size, and each token predicts one point-set. The first
point-set in the predicted sequence describes the current state of the
kinematic chain in the given RGB image. The remaining point-sets are
predictions of the future states of the kinematic chain in the next $T-1$
time-steps, i.e., future robot movements. The network is trained to match each
point-set to the ground-truth kinematic structure at the corresponding time-
step by minimizing the earth-moving distance. The same process is applied to
multiple RGB images taken from different viewpoints. The network is also
trained to make the future point-sets predicted from different viewpoints
consistent with each other. The kinematics tokens are concatenated with the
text tokens as query tokens. The concatenated query and vision tokens are
forwarded to our proposed multi-view dual attention block, to interweave the
information between query and vision tokens and across multiple camera
viewpoints.

The multi-view dual attention block performs dual attentions between query and
vision tokens within each single view, indicated by the green blocks in Figure
3, and multi-view attentions for query and vision tokens independently but
across multiple camera viewpoints, shown as the blue blocks in Figure 3. This
dual attention enables the query token to learn the visual kinematic chain
movements. The multi-view attention enables a spatially-consistent prediction
across multiple viewpoints. The output kinematics tokens are projected to $T$
point-sets using a small MLP network. Each point-set contains $N$ 2D points on
the input image. The VKT architecture is illustrated in Figure 3.

Compared to existing robot learning transformers [18, 17], our proposed VKT is
convolution-free and built solely with attention layers. Therefore, VKT
supports an arbitrary number of camera viewpoints, which is an important
advantage as camera setups are different among various environments.

Figure 4: Projecting the kinematics tokens into robot actions with 1D
convolution.

Deploying to a Specific Environment. We freeze the trained VKT, drop the
point-set prediction branch, and only use the kinematics tokens. We apply a 1D
convolution to project the $T$ kinematics tokens into low-level robot actions
for the next $T$ time-steps. Unlike the VKT backbone which is trained from
multiple environments, the 1D convolution head is trained separately for each
environment. The corresponding proposed architecture is illustrated in Figure
4.

Table 1: Comparison of Behavioral Cloning Transformer (BCT) and our Visual Kinematics Transformer (VKT) on Calvin and RLBench. BCT and VKT share the same architecture except that BCT directly predicts the robot’s end-effector movements. VKT outperforms BCT in isolated environments (specialized) and retains a similar performance when trained in both environments (general).  Environment | Metric | Specialized Agent | General Agent
---|---|---|---
BCT | VKT (Ours) | BCT | VKT (Ours)
Calvin [28] | Success Length $\uparrow$ | 1.1 | 1.173 | 0.334 | 1.083
RLBench [29] | Success Rate $\uparrow$ | 36.4% | 55.5% | 24.3% | 61.7%

##  4 Experiments

We evaluate VKT on Calvin [28] and RLBench [29], two standard benchmarks in
language-conditioned multi-task manipulation learning. Further, we evaluate
VKT in a subset of Open-X Embodiment [4] and real robot experiments. We aim to
answer the following questions: (1) Can visual kinematics forecasting improve
specialized agents in each environment? (2) Can VKT serve as a strong general
agent in multiple environments? (3) Can VKT be efficiently trained with real-
world demonstrations and solve manipulation tasks with real robots?

Setup. We adopt the setting of RT-X [5] for imitation learning from RGB images
without depth. RLBench contains various task categories, such as pick-and-
place, tool use, high-precision operations, screwing, tasks with visual
occlusion, and long-term manipulation. RLBench provides five RGBD cameras. We
use the RGB signal from the front, left, and right cameras. We choose 17 tasks
based on camera visibility and the task categorization of Hiveformer [30]. We
use 100 recorded trajectories per task for training, evaluate the trained
agent on 25 independent trials, and report the average success rate. Calvin
requires policies to accomplish multiple tasks sequentially in each episode.
We follow the convention of [28] of training on the pre-collected
demonstrations and reporting the average number of successful steps per
episode (success length). Examples of Calvin and RLBench tasks are shown in
Figure 5. Additional visualizations are provided in the supplementary
material.

Implementation Details. We use CLIP ViT-B/16 [27] as text and vision encoders.
We resize all RGB inputs to $224\times 224$. We use 8 multi-view dual
attention blocks ($L=8$) and predict visual kinematics chains for a horizon of
20 time-steps ($T=20$) and set size $D$ to $512$ in all experiments.

###  4.1 Learning from Multiple Environments

####  4.1.1 Simulation

Figure 5: Predicted visual kinematic chains returned by our VKT for different
robots in Calvin (left) and RLBench (right). The predicted kinematic chain of
the current frame is colored in red and the forecast one for the next time-
step is in blue. Videos are included in the supplementary material.

Baselines. We compare the performance of our visual kinematics transformer
(VKT) and a behavioral cloning transformer (BCT). The BCT shares the same
neural architecture with VKT. The difference is that BCT is directly trained
to predict 6-DoF end-effector poses and gripper states. In contrast, VKT is
trained to forecast visual kinematic chains first, then project the kinematics
tokens to robot actions using a convolution head with the backbone frozen, as
shown in Figure 4. The specialized agent is trained on each environment
separately. The general agent is trained using demonstrations from both
environments. We use the 6-DoF end-effector pose and binary gripper states as
the robot action space for both Calvin and RLBench. Note that RLBench requires
an additional gripper state for collision avoidance. We use two parallel
convolution heads for the BCT general agent to predict the low-level robot
actions for each environment.

Table 2: Performance on RLBench. We report the success rates on each task and
the average overall success rate. Success rates are measured from five
independent runs.

Method | Avg. | Put | Reach | Turn | Slide | Open | Money | Place | Sweep
---|---|---|---|---|---|---|---|---|---
Success | In Drawer | Drag | Tap | Blocks | Drawer | In Safe | Wine | To Pan
Specialized Agent | BCT | 36.4 | 1.6 | 77.6 | 63.2 | 12.0 | 20.8 | 44.8 | 44.0 | 5.6
VKT (Ours) | 55.5 | 52.0 | 88.8 | 40.8 | 22.4 | 56.0 | 76.8 | 31.2 | 48.0
General Agent | BCT | 24.3 | 0.0 | 28.4 | 68.8 | 8.0 | 8.6 | 8.0 | 4.0 | 0.0
VKT (Ours) | 61.7 | 48.8 | 77.6 | 45.6 | 12.8 | 73.6 | 81.6 | 89.6 | 50.4
|  | Meet | Phone | Lid | Close | Close | Ball | Push | Lift | Close
|  | Off Grill | On Base | Off Pan | Microwave | Box | In Hoop | Button | Block | Jar
Specialized Agent | BCT | 32.0 | 37.6 | 47.2 | 93.6 | 84.0 | 33.6 | 21.6 | 0.0 | 0.0
VKT (Ours) | 41.6 | 68.0 | 92.0 | 64.0 | 93.6 | 85.6 | 64.0 | 8.8 | 10.4
General Agent | BCT | 21.6 | 0.0 | 52.0 | 80.0 | 76.8 | 36.4 | 20.0 | 0.0 | 0.0
VKT (Ours) | 33.6 | 72.8 | 98.4 | 92.0 | 84.8 | 93.6 | 69.6 | 11.2 | 12.8

Results. Table 1 summarizes the comparison between VKT and BCT. By forecasting
the visual kinematics chains, VKT outperforms BCT in both specialized and
general agents. This indicates that the kinematics tokens capture sufficient
information about the robot’s action, despite not seeing any low-level robot
actions. When trained in both environments, the general BCT agent suffers a
catastrophic performance drop (1.1 $\rightarrow$ 0.33, $36.4\%\rightarrow
24.3\%$). In contrast, the general VKT agent only has a minor drop in Calvin
($1.173\rightarrow 1.083$) but its performance increases ($55.5\%\rightarrow
61.7\%$) in RLBench. This indicates that by learning actions in the image
planes, training becomes more robust to distribution shifts. This is further
verified by our real-robot experiments in Section 4.2. We list the per-task
and per-step success rates of RLBench and Calvin in Table 2 and 6. Overall,
VKT outperforms BCT except a lower success rate in the “turn tap” task. This
task requires some specific end-effector rotations which involve large
variations in the Euler angles of the end-effector but small changes in the
visual kinematic chain from all camera viewpoints. We discuss these
limitations and future directions in Section 5.

####  4.1.2 Open-X Embodiment

Setup. Open-X Embodiment (Open-X) [5] is a collection of video demonstrations
of real-world manipulation tasks crowd-sourced from diverse environments.
Since Open-X does not provide robot models such as URDF, and camera parameters
are not available in most datasets, we select four datasets from Open-X and
manually annotate the robot joint positions from $80$ demonstration videos
with $1,1359$ frames. Using the annotated videos, we train a single VKT for
the four datasets. Then, we train the network to predict low-level robot
actions as in Figure 4. However, instead of freezing the VKT, we continue
finetuning the entire network with the trained VKT as initial weights. For the
BC transformer, we use the same training setup but randomly initialize the
network.

Figure 6: Visual kinematic chains of different robots predicted by our VKT in
Open-X Embodiment (left) and a comparison of the relative precision of the
behavior cloning transformer (BCT) and our visual kinematics transformer (VKT)
on the four selected datasets.

Results. Figure 6 shows qualitative examples of the predicted kinematic chains
(left) and compares the relative precision of BCT and VKT on the validation
set (right). Let $\mathcal{L}_{vkt}$ and $\mathcal{L}_{bct}$ denote the L1
error of VKT and BCT, the relative precision of VKT is computed as
$|\mathcal{L}_{vkt}-\mathcal{L}_{bct}|/\mathcal{L}_{vkt}$. Figure 6 shows that
our VKT learns to predict kinematic chains from diverse real-world data, and
visual kinematics forecasting can serve as a pre-training objective to improve
imitation learning.

###  4.2 Real Robot Experiments

Figure 7: Visual chains predicted by VKT in a language-conditioned pick-and-
place task with a real robot. The chains of pick and place actions are in red
and blue, respectively.

Table 3: Success rates of VKT for solving the real-world language-conditioned
pick-and-place task.

Object | BCT | VKT (Ours)
---|---|---
Overall (%) | 41.5 | 69.2
Crackers | 3/7 | 4/7
Cup | 3/7 | 4/7
Ball | 3/8 | 7/8
Chips | 1/7 | 3/7
Brick | 5/8 | 6/8
Tomato Can | 4/7 | 6/7
Mustard | 0/7 | 4/7
Sugar Box | 4/7 | 6/7
Orange | 4/7 | 4/7

Setup. We evaluate our VKT in a real-world language-conditioned pick-and-place
task. The task input includes an RGB image and text instructions such as “put
the cup into the blue box”. We adopt nine YCB objects such as balls, chips,
and two boxes as place targets. Examples of this task and the forecasted
kinematic chains are shown in Figure 7. We implement a scripted policy using a
DINOv2-based object detector [31] to collect 90 demonstrations as training
data, with only 5 demonstrations for each object-container pair. We
automatically compute the ground-truth visual kinematic chain with the robot’s
URDF and the calibrated camera parameters. We use a Kuka LBR iiwa robot. For
the convolution head, we train the model to directly predict the pick and the
place poses instead of predicting the full trajectory for simplicity, and use
MoveIt [32] for path generation.

Results. Table 3 shows the success rates of VKT in the real-world task. The
overall success rate is the average of 65 independent runs with different
object layouts. The BC Transformer (BCT) shares the same training and
evaluation setup with VKT, but VKT outperforms BCT by a significant margin.
Further, we discover that visual kinematic forecasting enables the use of 2D
image augmentation for manipulation learning because the action space also
resides in image planes, which is important to prevent overfitting in our
experiment due to the small training set. In comparison, existing manipulation
learning augmentations are applied in the 3D space [33], which requires depth
that is less available than RGB data, and has fewer categories than 2D image
augmentations [34].

###  4.3 Ablation Studies

Table 4: General agent’s performance on forecasting the full kinematic chain
vs. only the end-effector

Configuration | Calvin | RLBench
---|---|---
Success Length $\uparrow$ | Success Rate $\uparrow$
BCT | 0.334 | 24
VKT (Ours) | Full Chain | 1.083 | 55.5
End-Effector | 1.06 | 54.1

Table 5: Effect of multiple viewpoints in RLBench

Configuration | Success Rate
---|---
Single View | 26.8%
Multi View | 55.5%

Table 6: Performance on Calvin. Calvin requires agents to accomplish 5 tasks sequentially at each episode. We measure the step-level success rate and the average number of successful steps (success length). Success rates are measured from 1000 episodes. Method | Success Length $\uparrow$ | Success rates at each step (%)
---|---|---
1 | 2 | 3 | 4 | 5
Specialized Agent | BCT | 1.1 | 59 | 31.7 | 12.7 | 5.2 | 1.4
VKT (Ours) | 1.173 | 60.5 | 32.5 | 14.7 | 7.1 | 2.5
General Agent | BCT | 0.334 | 26.8 | 5.3 | 1.2 | 0.1 | 0
VKT (Ours) | 1.083 | 58.3 | 29 | 13.5 | 5.6 | 1.9

Predicting the Full Chain versus Predicting Only the End-Effector. We compare
the effects of forecasting the full chain and only the end-effector in Table
5. Table 5 shows that forecasting the full chain improves the performance on
both Calvin and RLBench. However, despite the inability to capture rotations,
end-effector VKT also delivers stronger performance as a general agent than
BCT. This indicates that a key to building general agents across diverse
environments is to predict robot movements directly in the image plane.
Furthermore, the end-effector predictions can be useful in special cases when
the robot’s body is largely invisible due to occlusions, or when robot models
are unavailable. This suggests a broader applicability for visual kinematics
forecasting.

Multi-View versus Single-View Predictions. We study the impacts of multi-view
forecasting in Table 5. We remove the multi-view attention layers to create
the single-view VKT. We average the low-level robot actions predicted from
each view as the final prediction. Table 5 shows that VKT improves its 3D
understanding significantly with our multi-view dual attention block by
predicting consistent kinematic chains in different camera viewpoints.

##  5 Discussion and Conclusion

We have shown that visual kinematics forecasting of quasi-static robot
movements can be a building block technique for developing general agents
across diverse robot learning environments. However, there are still some open
questions that are worth discussing. For example, we did not utilize wrist
cameras in VKT because it does not take images of the robot. This can be
resolved with RGBD wrist camera and virtual view projection [18]. Also, some
public real-robot datasets lack robot URDF files and reliable camera
parameters, which creates several challenges for applying our method. Finally,
it is worth noting that one can re-purpose VKT to forecast not only robot
movements, but also human poses and movements in order to leverage abundant
human demonstration videos for robot learning.

## References

  * Bommasani et al. [2021] R. Bommasani, D. A. Hudson, E. Adeli, R. Altman, S. Arora, S. von Arx, M. S. Bernstein, J. Bohg, A. Bosselut, E. Brunskill, et al.  On the opportunities and risks of foundation models.  _arXiv preprint arXiv:2108.07258_ , 2021.
  * Shridhar et al. [2022] M. Shridhar, L. Manuelli, and D. Fox.  Cliport: What and where pathways for robotic manipulation.  In _Conference on robot learning_ , pages 894–906. PMLR, 2022.
  * Shridhar et al. [2023] M. Shridhar, L. Manuelli, and D. Fox.  Perceiver-actor: A multi-task transformer for robotic manipulation.  In _Conference on Robot Learning_ , pages 785–799. PMLR, 2023.
  * Brohan et al. [2023] A. Brohan, N. Brown, J. Carbajal, Y. Chebotar, X. Chen, K. Choromanski, T. Ding, D. Driess, A. Dubey, C. Finn, et al.  Rt-2: Vision-language-action models transfer web knowledge to robotic control.  _arXiv preprint arXiv:2307.15818_ , 2023.
  * Padalkar et al. [2023] A. Padalkar, A. Pooley, A. Jain, A. Bewley, A. Herzog, A. Irpan, A. Khazatsky, A. Rai, A. Singh, A. Brohan, et al.  Open x-embodiment: Robotic learning datasets and rt-x models.  2023\.
  * Flamary et al. [2021] R. Flamary, N. Courty, A. Gramfort, M. Z. Alaya, A. Boisbunon, S. Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, L. Gautheron, N. T. Gayraud, H. Janati, A. Rakotomamonjy, I. Redko, A. Rolet, A. Schutz, V. Seguy, D. J. Sutherland, R. Tavenard, A. Tong, and T. Vayer.  Pot: Python optimal transport.  _Journal of Machine Learning Research_ , 22(78):1–8, 2021.  URL http://jmlr.org/papers/v22/20-451.html.
  * Rubner et al. [2000] Y. Rubner, C. Tomasi, and L. J. Guibas.  The earth mover’s distance as a metric for image retrieval.  _International journal of computer vision_ , 40:99–121, 2000.
  * Feydy [2020] J. Feydy.  Geometric data analysis, beyond convolutions.  _Applied Mathematics_ , 2020.
  * Dasari et al. [2020] S. Dasari, F. Ebert, S. Tian, S. Nair, B. Bucher, K. Schmeckpeper, S. Singh, S. Levine, and C. Finn.  Robonet: Large-scale multi-robot learning.  In L. P. Kaelbling, D. Kragic, and K. Sugiura, editors, _Proceedings of the Conference on Robot Learning_ , volume 100 of _Proceedings of Machine Learning Research_ , pages 885–897. PMLR, 30 Oct–01 Nov 2020.  URL https://proceedings.mlr.press/v100/dasari20a.html.
  * Huang et al. [2020] W. Huang, I. Mordatch, and D. Pathak.  One policy to control them all: shared modular policies for agent-agnostic control.  In _Proceedings of the 37th International Conference on Machine Learning_ , pages 4455–4464, 2020.
  * Kurin et al. [2020] V. Kurin, M. Igl, T. Rocktäschel, W. Boehmer, and S. Whiteson.  My body is a cage: the role of morphology in graph-based incompatible control.  In _International Conference on Learning Representations_ , 2020.
  * Gupta et al. [2021] A. Gupta, L. Fan, S. Ganguli, and L. Fei-Fei.  Metamorph: Learning universal controllers with transformers.  In _International Conference on Learning Representations_ , 2021.
  * Todorov et al. [2012] E. Todorov, T. Erez, and Y. Tassa.  Mujoco: A physics engine for model-based control.  In _2012 IEEE/RSJ international conference on intelligent robots and systems_ , pages 5026–5033. IEEE, 2012.
  * Reed et al. [2022] S. Reed, K. Zolna, E. Parisotto, S. G. Colmenarejo, A. Novikov, G. Barth-Maron, M. Gimenez, Y. Sulsky, J. Kay, J. T. Springenberg, et al.  A generalist agent.  _arXiv preprint arXiv:2205.06175_ , 2022.
  * Team et al. [2023] O. M. Team, D. Ghosh, H. Walke, K. Pertsch, K. Black, O. Mees, S. Dasari, J. Hejna, C. Xu, J. Luo, et al.  Octo: An open-source generalist robot policy.  In _Proceedings of Robotics: Science and Systems_ , Delft, Netherlands, 2023.
  * Chi et al. [2023] C. Chi, S. Feng, Y. Du, Z. Xu, E. Cousineau, B. Burchfiel, and S. Song.  Diffusion policy: Visuomotor policy learning via action diffusion.  _arXiv preprint arXiv:2303.04137_ , 2023.
  * Brohan et al. [2022] A. Brohan, N. Brown, J. Carbajal, Y. Chebotar, J. Dabis, C. Finn, K. Gopalakrishnan, K. Hausman, A. Herzog, J. Hsu, et al.  Rt-1: Robotics transformer for real-world control at scale.  _arXiv preprint arXiv:2212.06817_ , 2022.
  * Goyal et al. [2023] A. Goyal, J. Xu, Y. Guo, V. Blukis, Y.-W. Chao, and D. Fox.  Rvt: Robotic view transformer for 3d object manipulation.  _arXiv preprint arXiv:2306.14896_ , 2023.
  * Xu et al. [2022] C. Xu, S. Yang, T. Galanti, B. Wu, X. Yue, B. Zhai, W. Zhan, P. Vajda, K. Keutzer, and M. Tomizuka.  Image2point: 3d point-cloud understanding with 2d image pretrained models.  In _European Conference on Computer Vision_ , pages 638–656. Springer, 2022.
  * Zhu et al. [2023] H. Zhu, H. Yang, X. Wu, D. Huang, S. Zhang, X. He, T. He, H. Zhao, C. Shen, Y. Qiao, et al.  Ponderv2: Pave the way for 3d foundataion model with a universal pre-training paradigm.  _arXiv preprint arXiv:2310.08586_ , 2023.
  * Mendonca et al. [2023] R. Mendonca, S. Bahl, and D. Pathak.  Structured world models from human videos.  _arXiv preprint arXiv:2308.10901_ , 2023.
  * Belkhale et al. [2024] S. Belkhale, T. Ding, T. Xiao, P. Sermanet, Q. Vuong, J. Tompson, Y. Chebotar, D. Dwibedi, and D. Sadigh.  Rt-h: Action hierarchies using language.  _arXiv preprint arXiv:2403.01823_ , 2024.
  * Yuan et al. [2024] C. Yuan, C. Wen, T. Zhang, and Y. Gao.  General flow as foundation affordance for scalable robot learning.  _arXiv preprint arXiv:2401.11439_ , 2024.
  * Zhang and Boularias [2024] X. Zhang and A. Boularias.  One-shot imitation learning with invariance matching for robotic manipulation.  In _Proceedings of Robotics: Science and Systems_ , Delft, Netherlands, 2024.
  * Gu et al. [2023] J. Gu, S. Kirmani, P. Wohlhart, Y. Lu, M. G. Arenas, K. Rao, W. Yu, C. Fu, K. Gopalakrishnan, Z. Xu, et al.  Rt-trajectory: Robotic task generalization via hindsight trajectory sketches.  _arXiv preprint arXiv:2311.01977_ , 2023.
  * Cuturi [2013] M. Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  _Advances in neural information processing systems_ , 26, 2013.
  * Radford et al. [2021] A. Radford, J. W. Kim, C. Hallacy, A. Ramesh, G. Goh, S. Agarwal, G. Sastry, A. Askell, P. Mishkin, J. Clark, et al.  Learning transferable visual models from natural language supervision.  In _International Conference on Machine Learning_ , pages 8748–8763. PMLR, 2021.
  * Mees et al. [2022] O. Mees, L. Hermann, E. Rosete-Beas, and W. Burgard.  Calvin: A benchmark for language-conditioned policy learning for long-horizon robot manipulation tasks.  _IEEE Robotics and Automation Letters (RA-L)_ , 7(3):7327–7334, 2022.
  * James et al. [2020] S. James, Z. Ma, D. R. Arrojo, and A. J. Davison.  Rlbench: The robot learning benchmark & learning environment.  _IEEE Robotics and Automation Letters_ , 5(2):3019–3026, 2020.
  * Guhur et al. [2023] P.-L. Guhur, S. Chen, R. G. Pinel, M. Tapaswi, I. Laptev, and C. Schmid.  Instruction-driven history-aware policies for robotic manipulations.  In _Conference on Robot Learning_ , pages 175–187. PMLR, 2023.
  * Oquab et al. [2023] M. Oquab, T. Darcet, T. Moutakanni, H. Vo, M. Szafraniec, V. Khalidov, P. Fernandez, D. Haziza, F. Massa, A. El-Nouby, et al.  Dinov2: Learning robust visual features without supervision.  _arXiv preprint arXiv:2304.07193_ , 2023.
  * Coleman et al. [2014] D. Coleman, I. Sucan, S. Chitta, and N. Correll.  Reducing the barrier to entry of complex robotic software: a moveit! case study.  _arXiv preprint arXiv:1404.3785_ , 2014.
  * Mitrano and Berenson [2022] P. Mitrano and D. Berenson.  Data augmentation for manipulation.  _arXiv preprint arXiv:2205.02886_ , 2022.
  * Xu et al. [2023] M. Xu, S. Yoon, A. Fuentes, and D. S. Park.  A comprehensive survey of image augmentation techniques for deep learning.  _Pattern Recognition_ , 137:109347, 2023.


# Excitable actin dynamics and amoeboid cell migration

Nicolas Ecker<EMAIL_ADDRESS>Department of Biochemistry, University
of Geneva, 1211 Geneva, Switzerland  Department of Theoretical Physics,
University of Geneva, 1211 Geneva, Switzerland  Karsten Kruse
<EMAIL_ADDRESS>Department of Biochemistry, University of Geneva, 1211
Geneva, Switzerland  Department of Theoretical Physics, University of Geneva,
1211 Geneva, Switzerland  NCCR Chemical Biology, University of Geneva, 1211
Geneva, Switzerland

###### Abstract

Amoeboid cell migration is characterized by frequent changes of the direction
of motion and resembles a persistent random walk on long time scales. Although
it is well known that cell migration is typically driven by the actin
cytoskeleton, the cause of this migratory behavior remains poorly understood.
We analyze the spontaneous dynamics of actin assembly due to nucleation
promoting factors, where actin filaments lead to an inactivation of the
nucleators. We show that this system exhibits excitable dynamics and can
spontaneously generate waves, which we analyse in detail. By using a phase-
field approach, we show that these waves can generate cellular random walks.
We explore how the characteristics of these persistent random walks depend on
the parameters governing the actin-nucleator dynamics. In particular, we find
that the effective diffusion constant and the persistence time depend strongly
on the speed of filament assembly and the rate of nucleator inactivation. Our
findings point to a deterministic origin of the random walk behavior and
suggest that cells could adapt their migration pattern by modifying the pool
of available actin.

##  I Introduction

The ability of cells to migrate is one of their most fascinating
characteristics. During mesenchymal migration, cells persistently polarize and
adhere to the substrate, which leads to persistent directional motion FW'09 ;
PRB'10 . In contrast, during amoeboid migration, cells frequently change their
polarization and hence their direction of motion. They also adhere less
strongly to the substrate than cells during mesenchymal migration. Amoeboid
migration can be observed for the soil amoeba Dictyostelium discoideum and for
immune cells, for example, dendritic cells. The random walk performed during
amoeboid migration is an important aspect of immune cells’ task to scan the
organism for pathogens. The origin of the random polarization changes during
amoeboid migration is largely unknown RAV'16  and it is not clear to what
extent cells can control the characteristics of their random walk.

Molecular noise is an obvious candidate for generating random migration ASB'18
; MBA'20 . The processes involved in generating migration are indeed subject
to noise due to the stochastic nature of molecular reactions. However, these
stochastic events take place on length and time scales that are small compared
to those characteristic of cellular random walks. It is not obvious how cells
could influence the strength of this noise and hence their migration behavior.
Fluctuating external cues could also generate random walks. Indeed, cells
respond to a multitude of external signals, notably, chemical or mechanical
gradients, and adapt their migration accordingly. Here, the cells have a
certain degree of control as they can tune the strength of their responses.
However, cellular random walks have been observed in the absence of external
cues C'01 ; SP'03 ; A'08 . Finally, there is the possibility that cells
generate internal polarization cues, which would give them the maximal
possible control over their behavior. In this context, spontaneous actin
polymerization waves have been proposed to provide such internal cues SEF'20 .

Actin is an important constituent of the cytoskeleton, which drives cell
migration. It assembles into linear filaments called F-actin, with two
structurally different ends. This structural polarity of actin filaments is
exploited by molecular motors that transform the chemical energy released
during hydrolysis of adenosine-triphosphate (ATP) into mechanical work. The
assembly and disassembly of F-actin is regulated by various cofactors. For
example, formins and the Arp2/3 complex nucleate new filaments. Actin
depolymerizing factor (ADF)/cofilin, on the other hand, can promote their
disassembly. Interestingly, there is evidence for feedback between the actin
cytoskeleton and the activity of these regulatory cofactors. For example,
nucleation promoting proteins have been reported to be less active in regions
of high F-actin density BvD'17 ; BAG'09 . Such a feedback can lead to
spontaneous actin polymerization waves V'00 ; V'02 ; RV'12 ; WRO'06 ;
doubrovinski+kruseEPL . Such waves are present during migration V'00 ; WRO'06
; AM'13 ; WB'13 , and theoretical analysis has shown that they can be
sufficient to cause cell motility Doubrovinski+Kruse ; Dreher+Kruse ; SEF'20 ;
AM'13 .

From a physical point view, spontaneous actin polymerization waves are akin to
waves in excitable media. Early indications of this connection were given in
V'00 ; V'02 ; RV'12 . Further support came from the observation that actin
polymerization waves exhibit a refractory period EK'05 ; WRO'06 . More
recently, the actin network/cytoskeleton of D. discoideum was shown to be
poised close to an oscillatory instability WB'13 . The dynamics of excitable
systems is exemplified by the FitzHugh-Nagumo system, which is a very much
simplified version of the Hodgkin-Huxley equations describing action
potentials traveling along the axons of nerve cells.

In this work, we analyze the description of actin polymerization waves
proposed in Ref. doubrovinski+kruseEPL . We clarify its connection to the
FitzHugh-Nagumo system and characterize the waves it generates. Furthermore,
we use a phase-field approach SRL'10 ; ZSA'12  to study the impact of actin
polymerization waves on cell migration. Here, the phase field is an auxiliary
field that distinguishes between the inside and outside of a cell. We analyze
in detail a recently introduced current for confining proteins to the cell
interior SEF'20 . Finally, we explore the relation between the system
parameters and the characteristics of the random walks generated by chaotic
polymerization waves.

##  II Actin dynamics

In this section, we present the description of the actin cytoskeleton
developed in Refs. doubrovinski+kruseEPL ; Dreher+Kruse ; SEF'20 . After
establishing the dynamic equations, we discuss their relation to the FitzHugh-
Nagumo model (FHN) and show that oscillations and waves emerge spontaneously
in our system. Finally, we characterize the waves shape, length and
propagation velocity.

###  II.1 The dynamic equations

Amoeboid cell migration is driven by the actin cytoskeleton, which is mostly
concentrated in the actin cortex, a layer beneath the plasma membrane. The
cortex thickness is a few hundred nanometers CDP'13 ; CEF'17 ; JG'19  and thus
much smaller than the lateral extension of a cell ($>$10 $\mu$m). In this
work, we aim at describing the actin cytoskeleton adjacent to the substrate
and thus use a two-dimensional geometry.

We use the continuum description of Refs. doubrovinski+kruseEPL ; Dreher+Kruse
; SEF'20  for the actin dynamics, where the actin density is captured by the
field $c$. The alignement of actin filaments can lead to (local) orientational
order in the system. This effect is captured by the orientational order
parameter $\mathbf{p}$, which is similar to the nematic order parameter of
liquid crystals. In the dynamic equations, all terms allowed by symmetry up to
linear order and up to first order in the derivatives are considered, such
that

$\displaystyle\partial_{t}c$
$\displaystyle=-v_{a}\mathbf{\nabla}\mathbf{p}-k_{d}c+\alpha n_{a}$ (1)
$\displaystyle\partial_{t}\mathbf{p}$
$\displaystyle=-v_{a}\mathbf{\nabla}c-k_{d}\mathbf{p}.$ (2)

Here, $v_{a}$ is the average polymerization speed and $k_{d}$ an effective
degradation rate, see Fig. 1. Note, that this description neglects flows of
the actin network M'13  that could, for example, be generated by molecular
motors. We also neglect a possible diffusion term that would account for
fluctuations in the actin dynamics. We have checked that our results are not
affected qualitatively for sufficiently small diffusion constants. Equations
(1) and (2) can also be obtained by coarse-graining a kinetic description
Dreher+Kruse .

Figure 1: Schematic representation of the actin dynamics captured by Eqs.
(1)-(4). Blue circles represent inactive nucleators. They are spontaneously
activated at rate $\omega_{0}$, a process that is often associated with
membrane binding. The activation rate is enhanced by already active
nucleators, represented by green circles, which is captured by the parameter
$\omega$. Active nucleators generate new actin filaments (red) at rate
$\alpha$. The latter grow at velocity $v_{a}$ and spontaneously disassemble at
rate $k_{d}$. Furthermore, actin filament attract factors that inactivate
nucleators. This complex process, which can involve several proteins, is
captured by the rate $\omega_{d}$.

The last term of Eq. (1) is a source term that describes nucleation of new
actin filaments. For the conditions present in cells, new actin filaments
hardly form spontaneously. Instead, specialized proteins assist in this
process. Examples are members of the formin family or the Arp2/3 complex.
These proteins can be in an active or an inactive state and their spatial
distribution in a cell can change with time. In this way they can contribute
essentially to orchestrating the organization of the actin cytoskeleton. We
introduce the densities $n_{i}$ and $n_{a}$ to describe these actin nucleation
promoting factors - ’nucleators’ for short -, where the indices refer to the
inactive and active forms, respectively. Active nucleators generate new actin
at a rate $\alpha$, hence the form of the last term in Eq. (1).

The dynamic equations for the fields $n_{a}$ and $n_{i}$ capture their
transport by diffusion and their activation and inactivation dynamics. On the
time scales that are relevant for the dynamics we study in the remainder of
this work, nucleator synthesis and degradation can be neglected. Consequently,
the dynamic equations should conserve the number of nucleating proteins,
$\int_{A}(n_{a}+n_{i})dA=An_{\mathrm{tot}}=const$, where $A$ is the cell area
adjacent to the substrate. We write

$\displaystyle\partial_{t}n_{a}$ $\displaystyle=D_{a}\Delta
n_{a}+\omega_{0}\left(1+\omega n_{a}^{2}\right)n_{i}-\omega_{d}cn_{a}$ (3)
$\displaystyle\partial_{t}n_{i}$ $\displaystyle=D_{i}\Delta
n_{i}-\omega_{0}\left(1+\omega n_{a}^{2}\right)n_{i}+\omega_{d}cn_{a}.$ (4)

The diffusion constants for active and inactive nucleators are $D_{a}$ and
$D_{i}$, respectively. Spontaneous activation of nucleators occurs at rate
$\omega_{0}$. There is some experimental evidence for a positive feedback of
nucleator activation GW'14 , such that active nucleators promote the
activation of further nucleators. We capture this effect by the parameter
$\omega$. Nucleator deactivation can occur spontaneously. Furthermore it has
been proposed that nucleator deactivation can be induced by factors that are
recruited by actin filaments WMK'07 ; BAG'09 ; BvD'17 ; GW'14 . We assume that
the latter dominates GW'14  and neglect spontaneous deactivation. Actin
induced deactivation is controlled by the parameter $\omega_{d}$.

To fully determine the dynamics of the fields $c$, $\mathbf{p}$, $n_{a}$, and
$n_{i}$, Eqs. (1)-(4) have to be complemented by boundary conditions. In this
section, we use periodic boundary conditions to study the intrinsic actin
dynamics. Later we will add the presence of the cell membrane through a phase
field, see Sect. III.

In the following we use a non-dimensionalized version of the dynamic
equations. We scale time by $\omega_{0}^{-1}$ and space by
$\sqrt{D_{i}/\omega_{0}}$. We use the same notation for the rescaled
parameters as in Eqs. (1)-(4), such that the non-dimensionalization
corresponds to setting $\omega_{0}=1$ and $D_{i}=1$. Unless noted otherwise,
we use in the following the parameter values given in Table 1.

Table 1: Nondimensional parameter values used in this work unless indicated otherwise. The length and time scales are chosen such that the ensuing dynamics is comparable to that of immature dendritic cells SEF'20 . Parameter | Meaning | Value
---|---|---
$D_{a}$ | Diffusion constant of active nucleators | $4\cdot 10^{-2}$
$v_{a}$ | Effective actin polymerization speed |  $0.1$ \- $0.6$
$k_{d}$ | Effective filament degradation rate | $176$
$\omega$ | Cooperative binding strength of nucleators | $6\cdot 10^{-3}$
$\omega_{d}$ | Detachment rate of active nucleators | $0.1-0.6$
$\alpha$ | Actin polymerization rate | $588$
$n_{tot}$ | Average total nucleator density | 700
$L$ | System length | 1.3
$N_{g}$ | Number of grid points per dimension | 256
$\omega_{0}^{-1}$ | Time scale |  $91.6$ s
$\sqrt{D_{i}/\omega_{0}}$ | Length scale | 63.5 $\mu$m
$D_{\Psi}$ | Phasefield relaxation / surface tension coefficient | $5\cdot 10^{-3}$
$\kappa$ | Phasefield timescale modifier | $118$
$\epsilon$ | Area conservation strength | $8$
$\beta$ | Actin-membrane interaction coefficient | $5.75\cdot 10^{-3}$
$A_{0}$ | Mean cell area | $~{}0.083$

###  II.2 Spatially homogeneous solutions

Consider the case of homogeneous protein distributions. The constraint on the
nucleator density thus is $n_{a}+n_{i}=n_{\mathrm{tot}}=const$, where
$n_{\mathrm{tot}}$ is the average total nucleator density. According to Eq.
(2), the polarization field is decoupled from the other fields and will tend
to zero, $\mathbf{p}\to 0$, for $t\to\infty$. The remaining dynamic equations
become

$\displaystyle\partial_{t}c$ $\displaystyle=-k_{d}c+\alpha n_{a}$ (5)
$\displaystyle\partial_{t}n_{a}$ $\displaystyle=\left(1+\omega
n_{a}^{2}\right)(n_{\mathrm{tot}}-n_{a})-\omega_{d}cn_{a},$ (6)

where we have used $n_{i}=n_{\mathrm{tot}}-n_{a}$.

Equations (5) and (6) are reminiscent of the FitzHugh-Nagumo (FHN) system
FH'61 ; NAY'62 . In its general form, the latter is given by R'81 :

$\displaystyle\frac{1}{\epsilon}\partial_{t}w$ $\displaystyle=v-aw$ (7)
$\displaystyle\partial_{t}v$ $\displaystyle=-w+I+f(v).$ (8)

Equation (7) describes generation of the ’carrier’ $w$ by the ’driver’ $v$ and
degradation of $w$ with rate $a$. Here, $\epsilon\ll 1$ is a small parameter,
such that the dynamics of $w$ occurs on longer time scales than the one of
$v$. The second equation captures inhibition of $v$ by $w$ and $I$ is an
external stimulus. Finally, $f(v)$ describes a feedback of $v$ on its own
production: in general, it promotes generation of $v$ for small values of $v$,
whereas it inhibits its production for larger values of $v$.

A typical specific choice of $f$ is $f(v)=v-\frac{v^{3}}{3}$. In that case,
the system essentially depends only on the parameter $a$ and the external
stimulus $I$, because variations in $\epsilon$ do not affect the dynamics
qualitatively as long as $\epsilon$ is small. Although the stimulus can depend
on time, for the time being, we consider the case of constant $I$. Information
about the asymptotic behavior can be obtained by analyzing the nullclines in
phase space, that is, the curves defined by the respective conditions
$\dot{v}=0$ and $\dot{w}=0$ in the $(v,w)$-plane. Intersections of the two
nullclines correspond to fixpoints of which there are either one or three. In
the latter case, the system is bistable as two fixpoints are stable against
small perturbations, whereas the third is unstable, see Fig. 2A.

Figure 2: Phase space diagrams for spatially homogenous dynamics. A-C) Phase
space for the FitzHugh-Nagumo equations (7) and (8) with $a=2$, $I=0$ (A),
$a=0.04$, $I=2$ (B), and $a=0.4$, $I=2$ (C). D-F) Phase space for the dynamic
equations (5) and (6) with $k_{d}=5$, $\alpha=50$ (D), $k_{d}=50$,
$\alpha=400$ (E), and $k_{d}=80$, $\alpha=400$ (F). Other parameters as in
Table 1. In each case, the nullclines are shown in red, the vector fields as
blue arrowheads and an example trajectory in black. For the FHN equations, the
diagrams show a bistable case (A), a limit cycle (B) and an excitable case
(C). For Eqs. (5) and (6) we present limit cycles (D, E) and an excitable case
(F). For these equations, there is no bistable case.

In the case that there is one fixpoint, it can be stable or unstable against
small perturbations. If it is unstable, the system exhibits a limit cycle and
asymptotically oscillates, see Fig. 2B. In the opposite case, the FHN system
can present excitable dynamics, that is, even though the fixpoint is stable
against small perturbations, sufficiently large perturbations induce an
’excursion’ in phase space, before returning to the fixpoint, see Fig. 2C.
This behavior can be observed, when the intersection of the two nullclines is
left to the minimum or right to the maximum of the $v$-nullcline. If the
intersection is between the two extrema, the system spontaneously oscillates,
see Fig. 2B.

The similarity between the actin-nucleator dynamics, Eqs. (5) and (6), and the
FHN system becomes evident when choosing $c=w$, $n_{a}=v$, $\epsilon=\alpha$,
$a=k_{d}/\alpha$, $I=n_{\mathrm{tot}}$, and $f(v)=-v+\omega Iv^{2}-\omega
v^{3}$. The two dynamical systems differ in that the term $-w$ of Eq. (8)
corresponds to $-\omega_{d}vw$ in Eq. (6). Lastly, in contrast to $v$ and $w$
in the FHN system, which can take any real value, we now have $w\geq 0$ and
$0\leq v\leq n_{\mathrm{tot}}$. Note that, in the FHN system, $I$ is an
external signal and can depend on time, while the corresponding term
$n_{\mathrm{tot}}$ in the actin-nucleator system is a constant.

From the comparison between the actin-nucleator dynamics and the FHN system,
we see that the actin-nucleator dynamics is driven by the nucleators, whereas
actin is the carrier providing negative feedback. This is in agreement with
experimental observations RV'12 ; EK'05 . The similarity between the two
systems suggests that the actin-nucleator dynamics can also show oscillations
as well as excitable behavior. This is indeed the case as we discuss now. We
consider the case, where $\alpha$ is not a small parameter.

Let us now take a closer look at the nullclines. Analogously to $\dot{w}=0$
for the FHN system, $\dot{c}=0$ yields a linear relation between $c$ and
$n_{a}$ and the $n_{a}$-nullcline exhibits the characteristic S-shape of
$\dot{v}=0$. The nullclines of our system intersect exactly once in the region
$c\geq 0$ and $n_{a}\geq 0$, such that there is only one fixpoint
$(c_{0},n_{a,0})$, independently of the parameter values. To see this, note
first that the $c$-nullcline is a straight line through the origin. Now
consider the function $c(n_{a})$ defined by the nullcline $\dot{n}_{a}$. If
there were parameter values for which three intersection points existed, then
there would be some tangent to $c(n_{a})$ with a negative y-intercept $c_{y}$.
However, for any value $n_{a}\geq 0$ the value $c_{y}$ is given by

$\displaystyle c_{y}$ $\displaystyle=\omega
n_{a}^{3}+2n_{\mathrm{tot}}-n_{a},$ (9)

which is always positive as the number of active nucleators is bounded from
above by the total number of nucleators, $n_{\mathrm{tot}}\geq n_{a}$, proving
the above statement.

If the fixpoint is unstable against small perturbations, the system exhibits
oscillations as mentioned above, see Fig. 2D, E. In case, $(c_{0},n_{a,0})$ is
stable, the system can amplify a finite perturbation, but will eventually
return to the fixpoint,see Fig. 2F. before performing a linear stability
analysis of the fixpoint, we first obtain a physical picture of the necessary
conditions for an instability based on the nullclines.

The fixpoint can only be unstable, when the $n_{a}$-nullcline $c(n_{a})$
exhibits two extrema for $n_{a}>0$. Explicitly, the nullcline is given by

$\displaystyle c(n_{a})$ $\displaystyle=\frac{n_{\mathrm{tot}}-\omega
n_{a}^{3}+\omega n_{\mathrm{tot}}n_{a}^{2}-n_{a}}{\omega_{d}n_{a}}.$ (10)

Consequently, $\lim\limits_{n_{a}\rightarrow\infty}c(n_{a})=-\infty$ and
$\lim\limits_{n_{a}\rightarrow 0^{+}}c(n_{a})=+\infty$. To determine whether
the $n_{a}$-nullcline is monotonously decreasing, we consider the positive
roots of the derivative $c^{\prime}=\partial c/\partial n_{a}$. They are
determined by

$\displaystyle 0$ $\displaystyle=-n_{\mathrm{tot}}-2\omega n_{a}^{3}+\omega
n_{\mathrm{tot}}n_{a}^{2}.$ (11)

This equation always has a negative real solution. Two positive roots can only
exist if the discriminant of the polynomial is negative. This leads to $\omega
n_{\mathrm{tot}}^{2}>27$. In that case, the two real roots take the form

$\displaystyle n_{a}^{\pm}$
$\displaystyle=\frac{n_{\mathrm{tot}}}{6}\left(1\pm
2\sin\left[\frac{\pi-\sin^{-1}\left(1-\frac{54}{\omega
n_{\mathrm{tot}}^{2}}\right)}{3}\right]\right).$ (12)

The value of $n_{a}^{+}$ is always positive and $n_{a}^{-}$ is always
negative, because the argument of the sine function takes values between
$\pi/6$ and $\pi/2$. The second positive root is

$\displaystyle n_{a}^{0}$
$\displaystyle=\frac{n_{\mathrm{tot}}}{6}\left(1+2\sin\left[\frac{\sin^{-1}\left(1-\frac{54}{\omega
n_{\mathrm{tot}}^{2}}\right)}{3}\right]\right).$ (13)

In conclusion, the fixpoint $(c_{0},n_{a,0})$ is unstable and the system
oscillates for $\omega n_{\mathrm{tot}}^{2}>27$ and if
$n_{a}^{0}<n_{a,0}<n_{a}^{+}$.

We now turn to a linear stability analysis of the fixpoint. For the dominating
growth exponent $s$ of the perturbation, we find

$\displaystyle s$
$\displaystyle=\frac{a-k_{d}+\sqrt{(a-k_{d})^{2}-4\alpha\omega_{d}n_{a,0}}}{2},$

where $a=-1-3\omega n_{a,0}^{2}+2\omega
n_{\mathrm{tot}}n_{a,0}-\omega_{d}c_{0}$ only depends on ${k}_{d}/\alpha$. By
increasing the nucleation rate $\alpha$ while keeping ${k}_{d}/\alpha=const$
the nullcline remains unaffected. For $k_{d}>a$ the real part of the
eigenvalue becomes negative, leading to a stationary state. Thus, $k_{d}<a$ is
the last condition for the presence of oscillations in our system. The
oscillation frequency $\omega_{F}$ close to the instability can be estimated
from the imaginary part of the growth exponent $s$ of a small perturbation
through $\omega_{F}=\Im(s)=\sqrt{\alpha\omega_{d}n_{a,0}}$.

###  II.3 Wave solutions

After having analyzed the dynamic equations (1)-(4) for spatially homogenous
fields, we now turn to the general case and study the system in a domain of
size $L^{2}$ with periodic boundary conditions in the $x$\- and $y$-direction.
Then, the system can generate a variety of spatially heterogeneous solutions,
including planar traveling waves and stationary patterns, see Fig. 3 and
Supplementary Movies 1,2. In the following we will determine the parameter
region, in which these patterns exist and characterize the shape of planar
waves.

Figure 3: Snapshots of solutions for the actin concentration $c$ to Eqs.
(1)-(4) in two dimensions with periodic boundary conditions. A, B) Travelling
planar waves for $D_{a}=0.04$, $\omega_{d}=0.28$, $v_{a}=0.2$ (A) and
$D_{a}=0.04$, $\omega_{d}=0.32$, $v_{a}=0.44$ (B). Green arrows indicate the
direction of motion. The disclinations in (B) might heal after very long
times. C, D) Stationary Turing patterns for $D_{a}=0.04$, $\omega_{d}=0.45$,
$v_{a}=6.0$ (C) and $D_{a}=0.21$, $\omega_{d}=0.42$, $v_{a}=9.5$ (D). For
different initial conditions a pure hexagonal pattern of blobs can appear. All
other parameters as in Table 1.

####  II.3.1 Linear stability analysis

We start our analysis by investigating the stability of the homogenous steady
state against small spatially heterogeneous perturbations. The homogenous
state is characterized by $c(x)=c_{0}=\alpha n_{a}/k_{d}$,
$\mathbf{p}(x)=\mathbf{p}_{0}=0$, and $n_{i,0}=n_{\mathrm{tot}}-n_{a,0}$ with

$\displaystyle(1+\omega n_{a,0}^{2})n_{i,0}-\omega_{d}c_{0}n_{a,0}$
$\displaystyle=0.$ (14)

As shown above there is only one positive solution $n_{a,0}\leq
n_{\mathrm{tot}}$ to this equation, such that there is a unique homogenous
stationary state.

Consider $c(x,y,t)=c_{0}+\delta c(x,y,t)$ and similarly for the fields
$\mathbf{p}$, $n_{a}$, and $n_{i}$. Linearizing the dynamic equations with
respect to the steady state and expressing the perturbations in terms of a
Fourier series, $\delta
c=\sum_{n,m=-\infty}^{\infty}\hat{c}_{nm}e^{-i(q_{x,n}x+q_{y,m}y)}$ and
similarly for $\delta\mathbf{p}$, $\delta n_{a}$, and $\delta n_{i}$ with
$q_{x,n}=2\pi n/L$ and $q_{y,m}=2\pi m/L$, leads to

$\displaystyle\frac{d}{dt}\hat{c}_{nm}$
$\displaystyle=-iv_{a}\left(q_{x,n}\hat{p}_{x,nm}+q_{y,m}\hat{p}_{y,nm}\right)$
$\displaystyle\quad-k_{d}\hat{c}_{nm}+\alpha\hat{n}_{a,nm}$ (15)
$\displaystyle\frac{d}{dt}\hat{p}_{x,nm}$
$\displaystyle=-iv_{a}q_{x,n}\hat{c}_{nm}-k_{d}\hat{p}_{x,nm}$ (16)
$\displaystyle\frac{d}{dt}\hat{p}_{y,nm}$
$\displaystyle=-iv_{a}q_{y,m}\hat{c}_{nm}-k_{d}\hat{p}_{y,nm}$ (17)
$\displaystyle\frac{d}{dt}\hat{n}_{a,nm}$
$\displaystyle=-D_{a}\left(q_{x,n}^{2}+q_{y,m}^{2}\right)\hat{n}_{a,nm}+\left(1+\omega
n_{a,0}^{2}\right)\hat{n}_{i,nm}$ $\displaystyle\quad+2\omega
n_{i,0}n_{a,0}\hat{n}_{a,nm}-\omega_{d}(c_{0}\hat{n}_{a,nm}+\hat{c}_{nm}n_{a,0})$
(18) $\displaystyle\frac{d}{dt}\hat{n}_{i,nm}$
$\displaystyle=-\left(q_{x,n}^{2}+q_{y,m}^{2}\right)\hat{n}_{i,nm}-\left(1+\omega
n_{a,0}^{2}\right)\hat{n}_{i,nm}$ $\displaystyle\quad-2\omega
n_{i,0}n_{a,0}\hat{n}_{a,nm}+\omega_{d}(c_{0}\hat{n}_{a,nm}+\hat{c}_{nm}n_{a,0}).$
(19)

The solutions to these equations are of the form
$\hat{c}\propto\mathrm{e}^{s_{nm}t}$ etc, where $s_{nm}$ are the growth
exponents of the modes $(n,m)$. If $s_{nm}>0$, then a heterogeneous steady
state emerges. If instead, $\Re(s_{nm})>0$ and $\Im(s_{nm})\neq 0$, then an
oscillatory state, that is, either a standing or a traveling wave, can be
expected.

Our numerical solutions indicate that all instabilities in our system are
super-critical such that there is no coexistence of different states that are
not linked by a symmetry transformation. Close to the instability, the
wavelength $\lambda_{0}$ of the unstable determines the wave length of the
emerging pattern. This remains true in a large region beyond the instability,
see Fig. 4. The wave length depends only weakly on the actin assembly velocity
$v_{a}$, Fig. 4A, B, and not on the nucleator inactivation rate $\omega_{d}$,
Fig. 4D, E. It increases with the diffusion constant $D_{a}$, Fig. 4C, and
decreases with the cooperativity parameter $\omega$, Fig. 4F.

Figure 4: Wavelength as a function of system parameters. Orange dots represent
values obtained from numerical solutions in two spatial dimensions with
periodic boundary conditions ($L=1.0$), blue lines are the results of a linear
stability analysis, see Sect. II.3.1. Parameter values are $\omega_{d}=0.44$
(A), $\omega_{d}=0.48$ (B), $v_{a}=0.46$ and $\omega_{d}=0.43$ (C),
$v_{a}=0.32$ (D), $v_{a}=0.48$ (E), and $v_{a}=0.46$ and $\omega_{d}=0.43$
(F). All other parameters as in Table 1.

In contrast to the wave length, we only get a poor estimate of the wave’s
propagation velocity from the linear stability analysis. In the following we
use a variational ansatz to determine the wave form and propagation velocity
of plane waves.

####  II.3.2 Wave form

We start by rewriting the dynamic equations (1)-(4). First of all, we combine
the equations for the actin density $c$ and the polarization $\mathbf{p}$ to
obtain one equation for the density. Furthermore, we exchange $n_{i}$ for
$N=n_{a}+n_{i}$. Finally, we consider solutions in a reference frame moving
with the wave velocity $v$. We will use periodic boundary conditions with
period $\Lambda$. We thus arrive at

$\displaystyle 0$
$\displaystyle=\frac{v^{2}-v_{a}^{2}}{\Lambda^{2}}\partial_{x}^{2}c+\left(\frac{v}{\Lambda}\partial_{x}+k_{d}\right)\left(k_{d}c-\alpha
n_{a}\right)$ (20) $\displaystyle-vN$
$\displaystyle=\frac{1}{\Lambda}\partial_{x}N-\frac{1-D_{a}}{\Lambda}\partial_{x}n_{a}-vn_{\mathrm{tot}}$
(21) $\displaystyle-\frac{v}{\Lambda}\partial_{x}n_{a}$
$\displaystyle=\frac{D_{a}}{\Lambda^{2}}\partial_{x}^{2}n_{a}+\left(1+\omega
n_{a}^{2}\right)(N-n_{a})-\omega_{d}cn_{a},$ (22)

where we have scaled space by $\Lambda$, such that the period is equal to $1$,
see App. A.

Equations (20) and (21) are linear and can be solved as soon as $n_{a}$ is
known, see App. A. To solve the nonlinear Eq. (22) we make the following
ansatz for a right-moving wave in the interval $[-1/2,1/2]$

$\displaystyle n_{a}(a_{1},a_{2},a_{3},a_{4},x)$
$\displaystyle=\frac{a_{1}}{2}e^{-a_{2}x}(1+\tanh[a_{3}x])(1-2x)^{a_{4}},$
(23)

where $a_{1}$ to $a_{4}$ are variational parameters. We constrain $a_{2}$ and
$a_{3}$ to vary in the intervals $[5,15]$ and $[30,50]$, respectively, whereas
$a_{4}$ can take on the values $2$, $3$, $4$; we do not impose any constraints
on $a_{1}$. Note that the test function (23) does not fulfill the periodic
boundary condition. However, since $a_{2},a_{3}\gg 1$,
$n_{a}(a_{1},a_{2},a_{3},a_{4},\pm 1/2)\approx 0$.

In our ansatz, the active-nucleator density $n_{a}$ increases according to the
exponential polynomial $x^{a_{4}}\mathrm{e}^{a_{2}x}$ at the front of the
wave. In this region actin is nucleated and increases correspondingly. The
trailing region of the wave is defined by a decrease of the active nucleator
density according to $1+\tanh(a_{3}x)$. This decrease results from a threshold
actin concentration beyond which nucleator inactivation occurs at a higher
rate than nucleator activation. Due to the large value of $a_{3}$, the
nucleator density drops sharply to zero and also the actin density decays
exponentially in the trailing region. The corresponding decay length is
$v/k_{d}$, see App. A.

After solving the linear equations (20) and (21), we calculate an error by
integrating the difference between the left and the right hand sides of (22)
over the whole period:

$\displaystyle Err(a_{1},a_{2},a_{3},a_{4},v)$
$\displaystyle=\int_{-\frac{1}{2}}^{\frac{1}{2}}\left|F(n_{a},c,N)\right|dx$
$\displaystyle F(n_{a},c,N)$
$\displaystyle=\frac{v}{\Lambda}\partial_{x}n_{a}+\frac{D_{a}}{\Lambda^{2}}\partial_{x}^{2}n_{a}-\omega_{d}cn_{a}+\left(1+\omega
n_{a}^{2}\right)(N-n_{a}).$ (24)

Minimizing the error yields values for the variational parameters $a_{1}$ to
$a_{4}$, $v$, and $\Lambda$.

In Figure 5A, we compare a solution obtained by the variational ansatz and by
numerically solving the dynamic equations (1)-(4). The agreement is very good
with the largest deviations being present at the front of the wave. Similarly,
the parameter dependence of the wave speed is reproduced well by our
variational ansatz, Fig. 5B, C. The wave speed is essentially independent of
the actin polymerization speed $v_{a}$ as long as $v_{a}\lesssim 1$, which is
consistent with our earlier remark that the wave dynamics is driven by the
nucleator activity rather than actin assembly. Furthermore, the wave speed
increases with the parameter $\omega_{d}$ describing nucleator inactivation by
actin. Indeed, as $\omega_{d}$ increases, nucleators are more rapidly
inactivated, such that they become available for activation at the wave front.

Figure 5: Shape and velocity of traveling waves in one dimension. A) Actin and
active nucleator concentrations $c$ (green) and $n_{a}$ (orange) for
$v_{a}=0.8$ and $\omega_{d}=0.35$. Dots are from a numerical solution, solid
lines are obtained from the variational ansatz Eq. (23) and the solution Eq.
(39). B, C) Wave speed as a function of the actin polymerization velocity
$v_{a}$ (B) and the nucleator inactivation parameter $\omega_{d}$ (C). All
other parameters as in Table 1.

####  II.3.3 Stationary patterns

In addition to planar traveling waves, the dynamic equations (1)-(4) can also
produce stationary patterns, see Fig. 3C, D. These Turing patterns appear if
$v_{a}\gtrsim 1$ and consist either of ’blobs’ of high or low active nucleator
densities or of labyrinthine stripes of high active nucleator density. These
structures can coexist in the same system. Since our focus in this work is on
actin waves, we refrain from discussing these states further.

##  III Cell motility from actin polymerization waves

Having analyzed the intrinsic actin dynamics, we now turn to a
characterization of cell migration patterns emerging form spontaneous actin
waves. We start by introducing a phase-field approach for describing the
cellular domain. It contains a novel current for confining the nucleators to
the cell interior. We then describe migration patterns and study the
dependence of their characteristics on the system parameters.

###  III.1 Phase-field dynamics

Similar to previous work on cell motility, we use a phase-field approach to
define the dynamic cell shape SRL'10 ; ZSA'12 . A phase field is an auxiliary
scalar field with values ranging between 0 and 1, which are called the pure
phases of the system. We treat values of 0 as being outside of the cell and
values of 1 as being inside. The phase-field dynamics is given by SRL'10 ;
ZSA'12

$\displaystyle\partial_{t}\Psi$
$\displaystyle=D_{\Psi}\Delta\Psi+\kappa\Psi(1-\Psi)(\Psi-\delta)-\beta\mathbf{p}\cdot\mathbf{\nabla}\Psi,$
(25) where $\displaystyle\delta$
$\displaystyle=\frac{1}{2}+\epsilon\left(\int_{A}\Psi
dA^{\prime}-A_{0}\right).$ (26)

The term proportional to $\kappa$ derives from a free energy with minima at
the pure phases. They are separated by an energy barrier at $\Psi=\delta$.
Conservation of the cell area/volume can be achieved by adjusting the value of
$\delta$ as described in Eq. (26): The actual cell area is given by
$\int_{A}\Psi dA^{\prime}$, it’s target area by $A_{0}$. If the cell is bigger
than $A_{0}$, then $\delta>0.5$ such that the overall cell area shrinks and
vice versa. For sufficiently large values of $\kappa$, the transition between
the two pure phases is sharp.

The transition region between the two pure phases determines the position of
the cell membrane. Specifically, we implicitly define the location of the cell
membrane by all positions $\mathbf{r}$ with $\Psi(\mathbf{r})=0.5$. The term
proportional to $D_{\Psi}$ accounts for interfacial tension between the two
pure phases and thus the surface tension of the membrane. For cells, surface
tension of the membrane dominates its bending energy ZSA'12 , which we
neglect. Finally, the term proportional to $\beta$ describes the interaction
strength between the phase field and the actin network. The interaction is
always directed along the polarization vector, such that the membrane can be
pushed outwards or pulled inwards ZSA'12 . In our case we do not observe
pulling to the inside.

The dynamics of the actin network and the nucleators is confined to the cell
interior by multiplying the dynamic equations (1)-(4) by $\Psi$. Conservation
of the nucleators is an important aspect of these dynamic equations. Simply
multiplying the corresponding transport term by $\Psi$ violates conservation
of the total nucleator amount and also leads to nucleators leaking out of the
cell interior Dreher+Kruse . Here, we choose a different option and instead
modify the nucleator current at the position of the membrane. For a particle
density $n$, we write

$\displaystyle\partial_{t}{n}$ $\displaystyle=D(\Psi\Delta n-n\Delta\Psi)$
(27)
$\displaystyle=\mathbf{\nabla}(D\Psi\mathbf{\nabla}n)-\mathbf{\nabla}(Dn\mathbf{\nabla}\Psi).$
(28)

This term evidently conserves the total particle number. It can be interpreted
as a combination of scaling the diffusion constant with $\Psi$ and introducing
an inwards flux proportional to $D$ at the membrane. This suggests that the
expression is efficient for keeping the nucleators inside the cell. This is
indeed the case as can be seen by solving for the stationary state of Eq.
(28), which is given by $n\propto\Psi$.

In this context, it is also instructive to look at the discretized version of
the right hand side of Eq. (28). Using the discretized Laplacian $\triangle
n_{j}\equiv(n_{j-1}-2n_{j}+n_{j+1})h^{-2}$, where $h$ is the discretization
length, we get in one dimension:

$\displaystyle D(\Psi_{j}\triangle n_{j}\ -n_{j}\triangle\Psi_{j})$
$\displaystyle=D\frac{n_{j+1}\Psi_{j}+n_{j-1}\Psi_{j}-n_{j}\Psi_{j+1}-n_{j}\Psi_{j-1}}{h^{2}}.$
(29)

From this expression it is evident that nucleators can hop only to a site $j$
inside the cell, i.e., with $\Psi_{j}>0$, see Fig. 6A.

Figure 6: Polymerization waves in presence of a phase field. A) Schematic
comparison of the discretized diffusion in absence (left) and presence (fight)
of a phase-field, see Eq. (29). B) Phase diagram of migration patterns as a
function of the actin growth velocity $v_{a}$ and nucleator inactivation
parameter $\omega_{d}$. C-E) Example trajectories with cell outlines drawn at
8 equidistant points in time for $v_{a}=0.34$ and $\omega_{d}=0.45$ (diffusive
migration, C), $v_{a}=0.22$ and $\omega_{d}=0.38$ (random walk with straight
segments, D), and $v_{a}=0.46$ and $\omega_{d}=0.43$ (random walk with curved
segments, E). Scale bars correspond to a length of 0.3. Other parameters as in
Table 1.

In presence of the phase field, the dynamic equations are

$\displaystyle\partial_{t}{c}$ $\displaystyle=\Psi(\alpha
n_{a}-v_{a}\mathbf{\nabla}\cdot\mathbf{p})-k_{d}c$ (30)
$\displaystyle\partial_{t}{\mathbf{p}}$
$\displaystyle=-v_{a}\Psi\mathbf{\nabla}c-k_{d}\mathbf{p}$ (31)
$\displaystyle\partial_{t}{n_{a}}$ $\displaystyle=D_{a}(\Psi\Delta
n_{a}-n_{a}\Delta\Psi)+\Psi((1+\omega n_{a}^{2})n_{i}-\omega_{d}cn_{a})$ (32)
$\displaystyle\partial_{t}{n_{i}}$ $\displaystyle=\Psi\Delta
n_{i}-n_{i}\Delta\Psi-\Psi((1+\omega n_{a}^{2})n_{i}-\omega_{d}cn_{a}).$ (33)

For actin, the diffusion current can be neglected as argued above, and thus
its dynamics is unaffected by the modified diffusion introduced in Eq. (28).
However, since the actin concentration is not a conserved quantity and rapidly
degraded in the absence of nucleators, we chose the degradation term to act
also outside the cell interior to get rid of any actin that might have left
the cell.

###  III.2 Actin-wave induced cell trajectories

In Figure 6B we show the phase diagram of the different dynamics patterns of
the phase field’s center
$\mathbf{r}_{c}=\int\mathbf{r}\Psi(\mathbf{r})d^{2}\mathbf{r}$ as a function
of the parameters $v_{a}$ and $\omega_{d}$. Five different dynamic states can
be distinguished. Below a critical value of $\omega_{d}$, waves do not emerge
in the system and the center settles into a stationary state. The critical
value of $\omega_{d}$ depends only weakly on $v_{a}$. There is a second
critical value, such that the center $\mathbf{r}_{c}$ is again stationary if
$\omega_{d}$ is larger than this critical value.

Close to the critical values of $\omega_{d}$, the actin-nucleator system forms
a spiral wave, see Movie 3. These spirals are symmetric and do not deform the
phase field. They spin around a fixed point, which coincides with the center
$\mathbf{r}_{c}$. Since the dynamic equations are isotropic, solutions with
clockwise or counter-clockwise rotations coexist. As the value of $\omega_{d}$
is, respectively, further increased or decreased, the spiral loses its
symmetry. In this case, the motion of the center $\mathbf{r}_{c}$ becomes
erratic and can be described as a random walk.

Three different types of random walks can be identified. First, the center
$\mathbf{r}_{c}$ can exhibit diffusive dynamics, see Fig. 6C and Movie 4.
Second, it can perform a random walk, where straight segments along which the
cell moves with constant velocity alternate with segments of diffusive motion,
see Fig. 6D and Movie 5. Also in the third type of random walk the cell center
changes between two states, namely, diffusive or curved motion, see Fig. 6E
and Movie 6. Along the curved segments, the radius curvature typically varies,
but there are special cases, for which the radius of curvature along the
curved segments is constant and the same for all segments. Note that for all
kinds of random walk trajectories, the direction of motion after a diffusive
segment is uniformly distributed. Similarly, the handedness of a curved
segment is uncorrelated with that of the preceding segment.

For the erratic trajectories, the actin-nucleator dynamics is chaotic. For the
persistent random walk, states in which axisymmetric waves emanate from a
center with a fixed position within the cellular domain. During the diffusive
states, we observe spiral wave chaos instead. In the states corresponding to
curved segments, the waves are not axisymmetric, which leads to ’protrusions’
of the membrane and a turning of the cell axis. In case of the diffusive
trajectories, the actin-nucleator dynamics exhibits spiral chaos. The
deterministic dynamic equations are thus able to replicate salient migration
features of searching cells SEF'20 .

###  III.3 Dependence of migration characteristics on parameter values

The random walks discussed above fall into the class of persistent random
walks. For a persistent random walk, the velocity of the walker has a finite
time autocorrelation, that is, its magnitude and direction persist for a
characteristic time $\tau$. Note that there are several realizations of a
persistent random walk. In a run-and-tumble process, the walker exhibits
periods during which it moves along straight lines with constant speed. These
periods are interrupted by events during which the walker essentially does not
move but changes its direction. Another possibility is that the direction of
motion and the speed varies constantly in a smooth way. Inbetween these
extremes, the segments of a run-and-tumble motion shows continuous changes of
the velocity. In all cases, the mean square displacement $\langle
r^{2}(t)\rangle$ is given by $\langle
r^{2}(t)\rangle=4Dt+2(v\tau)^{2}(e^{-t/\tau}-1)$. Here, $v$ is the mean
velocity of the persistent period and $D$ is the diffusion constant describing
the effective diffusive behavior on very long time scales. In the following we
study, how the effective parameters $\tau$, $v$, and $D$ depend on our system
parameters.

As shown in Figure 7A-C, the persistence time $\tau$, the speed $v$ and the
diffusion constant $D$ initially increase with $v_{a}$ and then decrease for
larger values of $v_{a}$. The non-monotonous behavior of these quantities is a
consequence of two competing effects. To see this, let us first recall that
the wave speed does not increase with increasing $v_{a}$, Fig. 5B. However,
the polarization of the actin network does increase in this case as can be
read of directly from Eq. (2). Consequently, the interaction between the actin
field and the membrane gets stronger and the membrane deformations are more
pronounced. At the same time, the pronounced membrane deformations feed back
on the actin waves, which are getting less regular. Thus, the cell
polarization is less efficient, such that the periods of persistent migration
are effectively shortened. At the same time, the migration speed decreases
during these periods. This is confirmed by the mean instantaneous speed of the
cell centers, which are very similar to the effective speed $v$, see Fig. 7B.

As a function of the parameter $\omega_{d}$, we observe a transition from a
persistent to a diffusive random walk. Below the transition, the parameters
$v$ and $D$ increase with $\omega_{d}$. In contrast, the value of $\tau$
depends non-monotonically on $\omega_{d}$; it first increases and then
decreases. Above a critical value of $\omega_{d}$, we find $\tau=0$. For these
values, the diffusion constant varies only slightly with $\omega_{d}$ and is
two orders of magnitude smaller than for the persistent random walks. The
dependence of $v$ on $\omega_{d}$ is linear for the persistent random walks.
Note that the values of $v$ obtained from fitting the mean square displacement
for $\tau\approx 0$ are not meaningful. The mean instantaneous velocity is
again very similar to $v$ for the persistent random walks. In the diffusive
regime, it still grows linearly with $\omega_{d}$. This is in line with the
wave velocity, which increases with $\omega_{d}$, see Fig. 5C.

Figure 7: Effective parameters of random walk trajectories. A-C) Diffusion
constant $D$ (A), speed $v$ (B), and persistence time $\tau$ (C) as a function
of the actin polymerization speed $v_{a}$. D-F) As (A-C), but as a function of
the nucleator inactivation parameter $\omega_{d}$. Values were measured by
fitting a persistent random walk model to the mean square displacement (MSD)
of the respective trajectories. In (B) and (E), also the mean speed measured
directly on the trajectories is shown (orange squares). Other parameters as in
Table 1.

##  IV Discussion

In this work, we have shown that a deterministic, self-organized system
describing the actin assembly dynamics in living cells is capable of
generating cellular random walks akin to amoeboid migration Dreher+Kruse ;
SEF'20 . We elucidated its relation to excitable systems by a comparison with
the FitzHugh-Nagumo system and characterized in detail spontaneously emerging
traveling waves. We recall that the wave propagation speed is independent of
the actin polymerization velocity $v_{a}$, such that the waves are driven by
the nucleator dynamics and not the actin dynamics.

By coupling the actin dynamics to a phase field, we studied the impact of the
spontaneous actin dynamics on cell migration. In this context, we introduced a
new expression for the nucleator current in presence of a phase field, such
that nucleators are confined to the cell interior. In other phase-field
studies of cell migration, conservation of particle numbers is typically not
an issue and all material leaving the cell interior is simply quickly degraded
ZSA'12 . If nucleators are not conserved, for example, by replacing the
concentration of inactive nucetaors $n_{i}$ by a constant, then the density of
active nucleators diverges and waves are absent from the system. In Ref.
Dreher+Kruse , nucleators that had leaked out of the system were reintroduced
into the cell by homogenously distributing them in the cell interior. In
contrast, the current $-D(\Psi\nabla n-n\nabla\Psi)$ used in this work acts
locally. All phases reported in Ref. Dreher+Kruse  are recovered and also the
topologies of the phase spaces are the same in both systems with one notable
exception: whereas in the present work erratic migration occurred for larger
values of $v_{a}$ and $\omega_{d}$ than for persistent migration, it was the
opposite in Ref. Dreher+Kruse .

By analyzing the mean-squared displacement of the simulated cells, we
characterized their persistent random walks in terms of a diffusion constant,
a persistence time, and the cell speed. We linked these effective parameters
to the actin-polymerization speed $v_{a}$ and the strength $\omega_{d}$ of the
negative feedback of actin on nucleator activity. It showed that these
parameters had a strong effect on the effective diffusion constant and the
persistence time, whereas the cell speed varied only by a factor of two. This
suggests that by changing the pool of available actin monomers, cells can
control important aspects of their random walks. This might allow notably
cells of the immune system patrolling an organism for pathogens to adapt their
behavior to the tissue they reside in.

A negative feedback of actin filaments on the nucleator activity is essential
for the emergence of spontaneous actin-polymerization waves. In cells,
indirect evidence has been found for this negative feedback, but it remains to
disentangle the molecular interactions involved. They likely involve the
action of small GTPases, which also take part in the signal transduction
pathways that couple external stimuli to the intracellular actin dynamics. In
future work it will be interesting to couple the actin-nucleator system to
such external signals and study the ensuing dynamics.

Furthermore, it will be interesting to study in future work collective cell
migration driven by spontaneous actin-polymerization waves. Previous phase-
field studies revealed how steric interactions between cells can lead to
collective migration LZA'15 ; NG'16  and how topographic surface structures
influence this behavior WAZ'19 . In the context of our work, one might expect
interesting synchronization phenomena between actin waves in different cells.

###### Acknowledgements.

We thank Carles Blanch-Mercader for helpful discussions and the Swiss National
Science Foundation (grant 205321-175996) for financial support.

##  Appendix A Wave profile

In this appendix, we determine the actin and nucleator densities for a wave
traveling at velocity $v$.

###  A.1 Actin density

The actin density $c$ and the polarisation field $p$ are given by Eqs.
(1)-(4), which in one spatial dimension and after non-dimensionalization read

$\displaystyle\partial_{t}c$ $\displaystyle=-v_{a}\partial_{x}p-k_{d}c+\alpha
n_{a}$ (34) $\displaystyle\partial_{t}p$
$\displaystyle=-v_{a}\partial_{x}c-k_{d}p.$ (35)

Deriving Eq. (34) with respect to time, we can eliminate the field $p$ and
obtain a linear equation for $c$ with an inhomogeneity proportional to
$n_{a}$:

$\displaystyle\partial_{t}^{2}c+2k_{d}\partial_{t}c+k_{d}^{2}c-v_{a}^{2}\partial_{x}^{2}c$
$\displaystyle=\alpha(k_{d}+\partial_{t})n_{a}.$ (36)

This is the equation for a wave with speed $v_{a}$, internal friction with
$2k_{d}$ and a driving proportional to $k_{d}^{2}$. The source of the wave
depends on $n_{a}$ and its time derivative. We will assume that the active
nucleators move as a solitary wave with velocity $v$, that is,
$n_{a}(x,t)=n(x-vt)$.

In the reference frame moving with the nucleation wave speed $v$ and
normalized by the wavelength $L$, Eq. (36) becomes

$\displaystyle\frac{{v}^{2}-{v}_{a}^{2}}{(k_{d}L)^{2}}\partial_{x}^{2}c-\frac{2v}{k_{d}L}\partial_{x}c+c$
$\displaystyle=\frac{\alpha}{k_{d}}(1-\frac{v}{k_{d}L}\partial_{x})n_{a}.$
(37)

The homogeneous solution $c_{h}(x)$ to this equation can be written as

$\displaystyle c_{h}(x)$
$\displaystyle=e^{{\lambda}x}\left[\left(\frac{v_{0}}{{\lambda}_{a}}-\frac{c_{0}{\lambda}}{{\lambda}_{a}}\right)\sinh({\lambda}_{a}x)+c_{0}\cosh({\lambda}_{a}x)\right],$
(38)

where ${\lambda}={k_{d}Lv}/({v}^{2}-{v}_{a}^{2})$ and analogously
${\lambda}_{a}={k_{d}Lv_{a}}/({v}^{2}-{v}_{a}^{2})$. In the above equation,
the amplitude of the homogeneous solution is fixed by the conditions
$c_{h}(0)=c_{0}$ and $c^{\prime}_{h}(0)=v_{0}$.

The solution to the in-homogeneous Eq. (37) with the source term
$S(x)=\frac{\alpha
k_{d}L^{2}}{{v}^{2}-{v}_{a}^{2}}(1-\frac{v}{k_{d}L}\partial_{x})n_{a}$ is
obtained by the method of variation of constants. We write $c_{0}=AS(x)$ and
$v_{0}=AS^{\prime}(x)$, where $A$ is the Wronskian of our system and arrive at
the full solution

$\displaystyle c(x)$ $\displaystyle=\frac{\alpha
Lv}{v^{2}-v_{a}^{2}}\int_{0}^{1}n_{a}(x+\xi)e^{-\lambda\xi}\left(\cosh(\lambda_{a}\xi)-\frac{v_{a}}{v}\sinh(\lambda_{a}\xi)\right)d\xi,$
(39)

where $\lambda=v/(v^{2}-v_{a}^{2})$ and analogously for $\lambda_{a}$.

The solution corresponds to a fraction of $\frac{v-v_{a}}{2v}$ of the scaled
nucleator density decaying on a lengthscale of
$L_{-}=\frac{1}{\lambda-\lambda_{a}}$ and a fraction of $\frac{v+v_{a}}{2v}$
decaying with $L_{+}=\frac{1}{\lambda+\lambda_{a}}$. The decaying part of the
actin wave can be fitted perfectly with the single parameter
$a\left(\frac{v-v_{a}}{2v}e^{-\lambda_{-}x}+\frac{v+v_{a}}{2v}e^{-\lambda_{+}x}\right)$.
Note that the nucleation rate $\alpha$ has no effect on the shape of the wave,
but only affects its amplitude.

The solution for the polarization field $p$ is obtained by solving Eq. (35)
for $p(x,t)\equiv p(x-vt)$.

###  A.2 Total nucleator density

We now rewrite the dynamic equations (1)-(4) for the active and inactive
nucleator concentrations $n_{a}$ and $n_{i}$ in terms of the total nucleator
concentration $N=n_{a}+n_{i}$ and $n_{a}$. In one spatial dimension and after
non-dimensionalization, we have

$\displaystyle\partial_{t}n_{a}$
$\displaystyle=D_{a}\partial_{x}^{2}n_{a}+\left(1+\omega
n_{a}^{2}\right)(N-n_{a})-\omega_{d}cn_{a}$ (40) $\displaystyle\partial_{t}N$
$\displaystyle=\partial_{x}^{2}N+(D_{a}-1)\partial_{x}^{2}n_{a}$ (41)

In the reference frame of the traveling wave, (41) becomes

$\displaystyle-\frac{v}{L}\partial_{x}N$
$\displaystyle=\frac{1}{L^{2}}\partial_{x}^{2}N+\frac{D_{a}-1}{L^{2}}\partial_{x}^{2}n_{a}.$
(42)

Integrating once and determining the integration constant by integrating once
more over the entire system, we arrive at a first order equation for the total
amount of nucleators,

$\displaystyle\partial_{x}N+vLN$
$\displaystyle=vLn_{\mathrm{tot}}+(1-D_{a})\partial_{x}n_{a}$ (43)

with $n_{\mathrm{tot}}$ being the average total nucleator density.

Equation (43) implies that with a homogeneous total nucleator concentration
$N=const=n_{\mathrm{tot}}$, gradients in $n_{a}$ also vanish. Thus, a
heterogeneity in the total nucleator concentrations is necessary to observe
waves and wave propagation requires nucleator transport.

Furthermore, $D_{a}$ is a measure for how far active nucleators can diffuse
around the bulk of the wave while bound before detaching, on a time scale
proportional to the wave period $\tau$, thus affecting the wave length.
$D_{a}$ needs to be sufficiently smaller than $D_{i}$ to create a length scale
difference large enough to enable the formation of the bulk of the wave and
maintain the imbalance in total nucleator concentration, otherwise the
constant distribution of proteins is the only solution (as the wave length
grows too large, or the imbalance shrinks too much to be supported).

The solution to Eq. (43) is given by

$\displaystyle N(x)$
$\displaystyle=n_{\mathrm{tot}}+(1-D_{a})\left[n_{a}(x)-vL\left(\frac{e^{-\frac{vL}{2}}}{2\sinh(\frac{vL}{2})}\int_{-\frac{1}{2}}^{\frac{1}{2}}n_{a}(\xi)e^{vL(\xi-x)}d\xi+\int_{-\frac{1}{2}}^{x}n_{a}(\xi)e^{vL(\xi-x)}d\xi\right)\right]$
(44)

From this equation we see that there are no waves, when $D_{a}=1(=D_{i})$.

###  A.3 Active nucleator density

Using the solutions for $c$, Eq. (39), and $N$, Eq. (44), we arrive at a
single equation for the distribution of the active nucleators in the reference
frame moving at the wave speed $v$:

$\displaystyle\frac{D_{a}}{L^{2}}\partial_{x}^{2}n_{a}(x)+\frac{v}{L}\partial_{x}n_{a}(x)=$
$\displaystyle\frac{\omega_{d}\alpha
Lv}{v^{2}-v_{a}^{2}}n_{a}(x)\int_{0}^{1}n_{a}(x+\xi)e^{-\lambda\xi}\left[\cosh(\lambda_{a}\xi)-\frac{v_{a}}{v}\sinh(\lambda_{a}\xi)\right]d\xi$
$\displaystyle-\left[1+\omega n_{a}(x)^{2}\right]n_{i}(x),$ (45) where
$\displaystyle n_{i}(x)=$ $\displaystyle
n_{\mathrm{tot}}-D_{a}n_{a}(x)-\frac{(1-D_{a})vL}{e^{vL}-1}\int_{0}^{1}n_{a}(x+\xi)e^{vL\xi}d\xi$
(46)

is the distribution of inactive nucleators. This non-linear integro-
differential equation can be solved using the variational ansatz of Sect.
II.3.

##  Appendix B Movie captions

Movie 1: Example of a traveling wave solution to Eqs. (1)-(4) in two
dimensions with periodic boundary conditions for $v_{a}=0.44$,
$\omega_{d}=0.32$. Other parameters as in Table 1. Disclinations can take very
long times to heal.

Movie 2: Example of a Turing pattern generated by Eqs. (1)-(4) in two
dimensions with periodic boundary conditions for $v_{a}=6.0$,
$\omega_{d}=0.45$, $D_{a}=0.04$. Other parameters as in Table 1.

Movie 3: Symmetric spiral wave solution of Eqs. (30)-(33) for $v_{a}=0.25$ and
$\omega_{d}=0.325$. Other parameters as in Table 1. Colors indicate the actin
concentration, red line corresponds to $\Psi=0.5$.

Movie 4: Asymmetric spiral wave solution of Eqs. (30)-(33) for $v_{a}=0.225$
and $\omega_{d}=0.35$, leading to diffusive motion. Other parameters as in
Table 1. Colors indicate the actin concentration, red line corresponds to
$\Psi=0.5$.

Movie 5: Wave solution of Eqs. (30)-(33) for $v_{a}=0.4$ and $\omega_{d}=0.4$,
leading to a dynamics of the phase field’s center, where straight segments
alternate with diffusive segments. Other parameters as in Table 1. Colors
indicate the actin concentration, red line corresponds to $\Psi=0.5$.

Movie 6: Wave solution of Eqs. (30)-(33) for $v_{a}=0.48$ and
$\omega_{d}=0.43$ leading to persistent random walk of the phase field’s
center, where curved segments alternate with diffusive segments. Other
parameters as in Table 1. Colors indicate the actin concentration, red line
corresponds to $\Psi=0.5$.

## References

  * (1) Friedl P, Wolf K: Plasticity of cell migration: a multiscale tuning model. J Cell Biol 2010; 188: 11–19.
  * (2) Pankovà K, Rösel D, Novotnỳ M, Bràbek J: The molecular mechanisms of transition between mesenchymal and amoeboid invasiveness in tumor cells. Cell Mol Life Sci 2010; 67: 63–71.
  * (3) Raynaud F, Ambühl M, Gabella C, Bornert A, Sbalzarini IF, Meister JJ, Verkhovsky AB: Minimal model for spontaneous cell polarization and edge activity in oscillating, rotating and migrating cells. Nat Phys 2016; 12: 367–373.
  * (4) Moreno E, Flemming S, Font F, Holschneider M, C Beta, S Alonso et al: Modeling cell crawling strategies with a bistable model: From amoeboid to fan-shaped cell motion. Physica D 2020; 412: 132591.
  * (5) Alonso S, Stange M, Beta C: Modeling random crawling, membrane deformation and intracellular polarity of motile amoeboid cells. PLoS ONE 2018; 13: e0201977.
  * (6) Chung CY, Funamoto S, A Firtel R: Signaling pathways controlling cell polarity and chemotaxis,Trends Biochem Sci 2001; 26: 557–566.
  * (7) Sohrmann M, Peter M: Polarizing without a C(l)ue, Trends Cell Biol 2003; 13: 526.
  * (8) Asano Y, Nagasaki A, Uyeda TQP: Correlated waves of actin filaments and PIP3 in Dictyostelium cells. Cell Motil Cytoskeleton 2008; 65: 923–934.
  * (9) Stankevicins L, Ecker N, Terriac E, Maiuri P, Schoppmeyer R, Vargas P, Lennon-Duménil AM, Piel M, Qu B et al: Deterministic actin waves as generators of cell polarization cues. Proc Natl Acad Sci USA 2020; 117: 826–835.
  * (10) Bretschneider T, Anderson K, Ecke M, Müller-Taubenberger A, Schroth-Diez B, Ishikawa-Ankerhold HC, Gerisch G: The Three-Dimensional Dynamics of Actin Waves, a Model of Cytoskeletal Self-Organization, Biophys J 2009; 96: 2888–2900.
  * (11) Bement WM, Leda M, Moe AM, Kita AM, Larson ME, Golding AE, Pfeuti C, Su KC, Miller AL, Goryachev AB et al: Activator-inhibitor coupling between Rho signalling and actin assembly makes the cell cortex an excitable medium Nat Cell Biol 2015; 17: 1471–1483.
  * (12) Vicker MG: Reaction-diffusion waves of actin filament polymerization/depolymerization in Dictyostelium pseudopodium extension and cell locomotion. Biophys Chem 2000; 84: 87–98.
  * (13) Vicker MG: Eukaryotic cell locomotion depends on the propagation of self-organized reaction-diffusion waves and oscillations of actin filament assembly. Exp Cell Res 2002; 275: 54–66.
  * (14) Weiner OD, Rentel MC, Ott A, Jedrychowski M, Yaffe MB, Gygi SP, Cantley LC, Bourne HR, Kirschner MW: Hem-1 complexes are essential for Rac activation, actin polymerization, and myosin regulation during neutrophil chemotaxis, PLoS Biol 2006; 4: 186–199.
  * (15) Doubrovinski K, Kruse K: Cytoskeletal waves in the absence of molecular motors. Europhys Lett 2002; 83: 18003.
  * (16) Ryan GL, Petroccia HM, Watanabe N, Vavylonis D: Excitable actin dynamics in lamellipodial protrusion and retraction. Biophys J 2012; 102: 1493–1502.
  * (17) Allard J, Mogilner A: Traveling waves in actin dynamics and cell motility. Curr Opin Cell Biol 2013; 25: 107–115.
  * (18) Westendorf C, Negrete Jr. J, Bae AJ, Sandmann R, Bodenschatz E, Beta C: Actin system close to the onset of oscillations. Proc Natl Acad Sci USA 2013; 110: 3853–3858.
  * (19) Doubrovinski K, Kruse K, Cell motility resulting from spontaneous polymerization waves. Phys Rev Lett 2011; 107: 258103.
  * (20) Dreher A, Aranson IS, Kruse K, Spiral actin-polymerization waves can generate amoeboidal cell crawling. New J Phys 2014; 16: 18.
  * (21) Edelstein-Keshet L: Mathematical Models in Biology. Philadelphia: SIAM; 2005.
  * (22) Shao D, Rappel WJ, Levine H:Computational Model for Cell Morphodynamics. Phys Rev Lett 2010; 105: 108104.
  * (23) Ziebert F, Swaminathan S, Aranson IS: Model for self-polarization and motility of keratocyte fragments. J R Soc Interface 2012; 9: 1084.
  * (24) Clark AG, Dierkes K, Paluch EK: Monitoring actin cortex thickness in live cells. Biophys J 2013; 105: 570–580.
  * (25) Clausen MP, Colin-York H, Schneider F, Eggeling C, Fritzsche M: Dissecting the actin cortex density and membrane-cortex distance in living cells by super-resolution microscopy. J Phys D 2017; 50: 064002.
  * (26) Jasnin M, Beck F, Ecke M, Fukuda Y, Martinez-Sanchez A, Baumeister W, Gerisch G: The Architecture of Traveling Actin Waves Revealed by Cryo-Electron Tomography. Structure 2019; 27: 1211–1223.
  * (27) Marchetti MC, Joanny JF, Ramaswamy S, Liverpool TB, Prost J, Rao M, Simha RA: Hydrodynamics of soft active matter. Rev Mod Phys 2013; 85: 1143–1189.
  * (28) Graziano BR, Weiner OD: Self-organization of protrusions and polarity during eukaryotic chemotaxis. Curr Opin Cell Biol 2014; 30: 60–67.
  * (29) Weiner OD, Marganski WA, Wu LF, Altschuler SJ, Kirschner MW: An actin-based wave generator organizes cell motility. PLoS Biol 2007; 5: e221.
  * (30) FitzHugh R: Impulses and physiological states in theoretical models of nerve membrane. Biophys J 1961; 1: 445–466.
  * (31) Nagumo J, Arimoto S, Yoshizawa S: An active pulse transmission line simulating nerve axon. Proc IRE 1962; 50: 2061–2070.
  * (32) Rinzel J: Models in neurobiology. In: Enns RH, editor. Nonlinear Phenomena in Physics and Biology. New York: Plenum Press; 1981. pp. 345–367.
  * (33) Löber J, Ziebert F, Aranson IS: Collisions of deformable cells lead to collective migration. Sci Rep 2015; 5: 9172.
  * (34) Najem S, Grant M:Phase-field model for collective cell migration. Phys Rev E 2016; 93: 052405.
  * (35) Winkler B, Aranson IS, Ziebert F: Confinement and substrate topography control cell migration in a 3D computational model. Commun Phys 2019; 2: 82.