Confinement and substrate topography control cell migration in a 3D computational model

Cell movement in vivo is typically characterized by strong confinement and heterogeneous, three-dimensional environments. Such external constraints on cell motility are known to play important roles in many vital processes e.g. during development, differentiation, and the immune response, as well as in pathologies like cancer metastasis. Here we develop a physics-driven three-dimensional computational modeling framework that describes lamellipodium-based motion of cells in arbitrarily shaped and topographically structured surroundings. We use it to investigate the primary in vitro model scenarios currently studied experimentally: motion in vertical confinement, confinement in microchannels, as well as motion on fibers and on imposed modulations of surface topography. We find that confinement, substrate curvature and topography modulate the cell’s speed, shape and actin organization and can induce changes in the direction of motion along axes defined by the constraints. Our model serves as a benchmark to systematically explore lamellipodium-based motility and its interaction with the environment. During development, wound healing, differentiation or cancer metastasis, cells move continuously in heterogeneous environments, hence understanding how cell migration is controlled by confinement and by different substrate shapes is crucial to begin to build a first conceptual framework for cell motility. Here the authors develop a computational approach to systematically investigate the effect that a complex environment has on cell motion and speed.

T he motility of living cells and tissue is an equally important topic for biology and non-equilibrium physics. Individual and collective cellular migration are examples of a selforganized non-equilibrium state characterized by the transduction of chemical energy delivered by the metabolism into directed mechanical motion. Recent studies demonstrated that several aspects of the seemingly complex dynamics of migrating cells and tissues can be rationalized in the framework of relatively simple physics-driven models 1-4 . Due to the relative ease in observation, data analysis, and modeling, quantitative physics-driven studies of cell motion have in the past focused almost exclusively on cells crawling along flat substrates 5 . This has been often criticized by the biology community 6 , since in vivo cells never live in such an idealized world. In contrast, cells migrate in rather complex three-dimensional (3D) environments that in addition are strongly inhomogeneous and also reorganized by both forces and chemical action from the cells. For instance, both leukocytes 7 and metastatic cancer cells 8 move along or within blood vessels until they reach a target site, where they need to squeeze through epithelial layers, or extracellular matrix (ECM). These examples already involve a plethora of topographical and confining features, namely motion on curved substrates, in tubular and planar confinement, on fibers or fiber bundles, as well as through random 3D fiber networks.
A realistic in vivo situation, e.g., the movement of cancer cells inside living tumors, is prohibitively difficult to investigate, although progress has been achieved using novel microscopy techniques 9 . Consequently, it is important to study well-defined scenarios in a systematic manner, with the goal to deduce general guiding principles. In the last years, experimental studies on cells in well-defined artificial or reconstituted 3D migration assays have gained pace, see Paul et al. 10 for a recent review and Maiuri et al. 11 for a comparison of different confining situations. A large part of these works focuses on cells in reconstituted networks of collagen or synthetic hydrogels [12][13][14] . In addition, well-defined, microfabricated systems have been developed and employed to study confined cell motion, like microchannels 15,16 and vertical confinement in slit-like geometries 17,18 . Substrates with topographical features have a longer tradition in tissue engineering 19,20 , and are now also employed to study motility [21][22][23][24][25][26] , as is the motion on micro-or nano-scale fibers [27][28][29] . Although a direct mapping from these model systems to the full in vivo case is in no way direct or obvious, these systems are in well-controlled conditions that allow to isolate different generic features of topography and/or confinement and their influence on cell motion.
On planar substrates, most cells move by forming lamellipodia, thin sheets containing actin filaments that polymerize and branch towards the cell membrane, thereby exerting protrusion forces 30,31 . A long-standing question had been whether these motile structures dominating 2D motion also exist and are relevant in 3D. There is growing evidence that this is indeed the case for fibroblasts 32 and neutrophils 33 , at least under certain conditions. Therefore, we here develop a physics-based modeling framework that transfers the knowledge of 2D lamellipodia-based cell motion 2,34,35 to 3D settings. The model is kept as simple as possible. It nevertheless is able to reproduce a plethora of the phenomena of single cell movement in confinement and on curved or topographically structured substrates found in the literature and puts them into a common framework. In fact, compared to earlier 3D models that exclusively studied flat substrates 1,[36][37][38] (with the exception of ref. 39 on fiber networks) the elegant use of the phase field approach 40at the same time for the cell and for the surroundings it interacts withallows to model motion in, in principle, arbitrarily shaped and topologically non-trivial environments.

Results
Computational model. We first briefly describe the modeling framework. The developed model is a non-trivial 3D extension of the 2D phase field approach presented in ref. 2 , known to describe  the subcritical onset of cellular motion 41 and realistic cell shapes,  with model extensions being able to describe guided motion on  inhomogeneous substrates 42,43 and lamellipodial waves 44 . The approach is conceptually simple as it has only two dynamic variables: first, the cell's interface is described by the transition region of a dynamic 3D scalar phase field ρ(r, t) ϵ [0, 1] for a closed domain representing the cell. The phase field value ρ = 1 corresponds to the cell's interior, ρ = 0 to the outside, and the ρ = 1/2 isosurface is identified with the cell membrane. And second, inside this domain, a 3D vector field p(r, t) describes the mean orientation of actin filaments and the degree of their parallel ordering (direction and absolute value of p, respectively). Arbitrary substrates and confinement scenarios can be modeled by specifying two additional, static phase fields: Φ(r) describing steric exclusion and Ψ(r) restricting the actin generation to regions close to the substrate. They have the same overall shape but the characteristic decay length of the latter is chosen to be larger, to allow for substantial actin inside the cell.
The cellular motility machinery is implemented via its basic physical features in the dynamics of ρ(r, t) and p(r, t): actin is generated close to the membrane (i.e., the phase field interface) and induces a protrusion force via the polymerization ratchet mechanism 45 , leading to an effective interface advection rate. Together with the cellular adhesion to the substrate, and assuming that actin polymerization is localized close to the substrate and predominantly tangential to the local substrate plane, the cell then spreads and flattens. Myosin motor-induced contraction then induces a symmetry breaking 46 and overall cell polarity 41 , with a region of high p at the front identified as the lamellipodium, and the cell moves with constant speed, see Supplementary Movie 1 for a typical polarization event.
Specifically, the 3D model equations read In Eq. (1), D ρ determines the width of the diffuse interface (and the cell's surface tension 47 , see also ref. 48 for modeling more realistic membrane features). The second term implements relaxational dynamics in the phase field potential Since the order parameter dynamics is non-conserved (see e.g. the discussion in ref. 40 ), it needs to be supplemented by a volume conservation. We chose a simple implementation via a global constraint, Here 1 2 is the stationary point, μ the constraint's stiffness and the term in brackets is the difference between the cell's current volume and the prescribed (initial) volume V 0 . The term α p⋅ ∇ρ accounts for the pushing force exerted by actin ratcheting against the cell's boundary, with α an effective velocity (i.e., force over a friction coefficient with the substrate). Adhesion is implemented on a coarse-grained level 49,50 via the term κ∇Φ ⋅ ∇ρ with an adhesion strength parameter κ, but explicit adhesive dynamics can in principle be introduced 42,51 . The last term describes excluded volume interaction, with a strength λ, via the potential More information on the last two contributions can be found in refs. 47,50 .
The polarization dynamics has three contributions from actin turnover, see the first line of Eq. (2): a diffusion (or elastic) term, a source term proportional to the parameter β (related to the actin polymerization velocity) and discussed in more detail below, and a sink term describing degradation of actin, e.g., via depolymerization, with time scale τ being of the order of the actin turnover time (typically seconds). For the cell to be able to polarize and move, the front-back symmetry has to be broken. Two mechanisms, based on key experimental observations 46 and associated to the action of myosin motors in the actin network, were already implemented in ref. 2 : a term σ|p| 2 added to the parameter δ in Eq. (1), that accounts for acto-myosin contraction with contraction rate σ and that can be motivated by active gel theory 52 . And a second term, γ(∇ρ ⋅ p)p, that models the presence of an anti-parallel acto-myosin bundle at the cell's rear 2,41 and explicitly breaks the front-back symmetry. Finally, the last term in Eq. (2) explicitly suppresses p inside the substrate.
Let us now discuss the source term of actin polarization, i.e., the second term in Eq. (2). In lamellipodia, actin nucleation is known to occur predominantly at the cell membrane, where regulators such as Wiskott-Aldrich Syndrome protein (WASP) and Arp2/3 are located 30,31 . This effect was modeled in ref. 2 via localizing the actin source term at the cell's interface by using the phase field gradient ∇ρ. However, an important issue arises in modeling 3D cell movement: while the projected 2D area of spreading and moving cells has been studied extensively 53 , and has been explained by the prevalent forces 54 (involving protrusion, adhesion, traction, and membrane bending), it is to date unknown what determines the thickness of the lamellipodium, which typically always lies in the range 100-400 nm. Although there exist some hypotheses, e.g., involving curvature-sensitive membrane-embedded regulators of actin nucleation 55 , the mechanisms remain speculative.
On the level of phenomenological description, we hence employed the following strategy: actin polymerization is modeled anisotropically with respect to the local tangential plane defined by the substrate. We introduce the (tensorial) projection operator onto the local tangential plane,P ¼Î À b n b n withÎ the identity and b n ¼ ∇Φ=j∇Φj the normal vector to the substrate.Pb n ¼ 0 holds and hence the term β½ð1 À θÞPð∇ρÞ þ θ∇ρ in Eq. (1) implements a polymerization rate of β in the tangential plane and of βθ normal to the substrate. The parameter θ ϵ [0, 1] hence governs the degree of this anisotropy, within the two limits of θ = 0 for polymerization exclusively in the tangential plane and θ = 1 for the isotropic case. The second main addition compared to earlier 2D models is the assumption that the actin polymerization rate decays with the distance from the substrate. This is implemented via the phase field Ψ(r), and enters the actin source term as a common factor. It is motivated by the fact that actin polymerization relies on regulatory processes close to the membrane, and has already been employed in a similar form in ref. 1 .
Let us briefly discuss the relevant scales and typical parameters. Eqs. (1) and (2) are already written in rescaled units, as discussed in ref. 2 . We use 1 μm and 10 s as the typical length and time scale, since then all activity-related parameters are of order one. Namely, the typical actin depolymerization rate of ≈1s −1 and maximum actin polymerization velocity ≈0.1-0.3 μms −1 , i.e.,~20-60 monomers per second, amount to the used values of 1/τ 1 = 0.1 and β = 1-3. We chose α = β/3 since the cell's interface movement should be smaller then the bare polymerization velocitynot all polymerization events create propulsion, and there is friction. Typical myosin motor speeds are of order 0.1-0.2 μms −1 , implying σ and γ of order one. Note that the maximum center-of-mass (c.o.m.) velocities of the modeled cellsbeing a model outcomeare of order one and hence correspond to the typical maximum speed of keratocytes (0.1 μms −1 ) 45  In the following we mainly discuss two distinct cell types: first, ideal "lamellipodial-type" cells, whose polymerization kinetics lies perfectly in the tangential plane defined locally by the substrate (θ = 0). This type corresponds best to the well-studied keratocytes, for which it has been shown that the actin network in the lamellipodium is very close to two-dimensional 56 . With sufficiently high adhesion strength to ensure the cell's contact with the substrate, both stationary and motile states can be obtained when varying the remaining actin-(α, β) and myosinrelated (γ, σ) parameters. A cell of this type moving steadily on a planar susbstrate is exemplified in Fig. 1a. The second cell type is refered to as of (partial) "wall pushing-type", having finite θ. For this type, actin generated perpendicularly to the substrate is nonnegligible, leading to stronger contact of the cell with the substrate, even in the absence of strong adhesion. Note that a somewhat similar motility mode was considered in ref. 57 . Again the system is bistable with both stationary and motile states. We however found that for the same acto-myosin parameters wallpushing cells are less motile.
Inspired by the plethora of recent experiments 15-18,21-29 , we applied our 3D model to study cell migration in confinement (slits and channels), its sensitivity to substrate curvature including motion on fibers, as well as motion on planar substrates with topographical features, see Fig. 1b-f for a summary.
Motion in vertical confinement. We first considered the generic scenario of vertical confinement, i.e., a slit geometry, where the cell is sandwiched between two co-planar walls 17 and can accommodate its shape and migration mode only in the quasi-2D plane. For the upper wall, we investigated two different types: either it has properties identical to the lower substrate, i.e., the cell can adhere, spread and form lamellipodia ("adhesive wall"). Alternatively, the upper wall is inert and just presents a steric obstacle for the cell to which it can neither adhere nor form a lamellipodium ("passive wall"). The latter can be easily modeled by putting κ = 0 there and by attributing to the wall only a steric field Φ and no actin regulatory field Ψ.
We investigated confined cell motion using the following generic preparation protocol: a top wall (either active or passive) was approached from above to a cell moving on the lower substrate gradually, in a step-by-step fashion, letting the moving cell accommodate to the new environment before approaching further. To test whether the effect of the confinement is reversible, after a strong confinement state had been reached, the upper wall was retracted in a similar step-by-step fashion. This protocol is a generalization of the classical cell compression experiment 58 to moving cells and could be experimentally realized employing micromanipulation techniques, e.g. involving a magnetic upper wall.
The most interesting scenario appears for strongly driven cells, as exemplified for the lamellipodium-type cell in a symmetric slit (i.e., with adhesive top wall) by the snapshots shown in Fig. 2a corresponding to the red curves in Fig. 2b, see also Supplementary Movie 3. In this case, there is a window of gap widths exhibiting a bistability of motile modes: for the same parameters, the cell can attach to and form lamellipodia on both walls equally, see the blue cell in Fig. 2a, with the cell speed increasing with decreasing the slit width. Alternatively, the cell can attach only to one of the walls (the lower one, as the top wall was approached), see the green cell in Fig. 2a and moves with a speed that is independent of gap width and higher than for a cell attached to both walls, due to a larger lamellipodium. Upon further decreasing the gap width, finally the cell is forced to attach to the top wall when the latter is sufficiently close. Consequently, vertical confinement constitutes a geometry-induced motility hysteresis loop: the motion of a cell in a spatially inhomogeneous vertical confinement will strongly depend on the history of its path. In contrast, if the wall is passive, this behavior is absent (orange curve in Fig. 2b). Only for strong confinement the cell speeds up a little as its lamellipodium, i.e., driven area, slightly increases.
Wall-pushing cells show a similar overall behavior in vertical confinement as lamellipodium-type cells. For the same driving parameters as the respective lamellipodium-type cells, they are however effectively less driven and hence less motile. Instead of exhibiting hysteresis, sufficiently weakly driven cells rather get stopped when starting to feel the upper wall, see the blue curve in Fig. 2b. Upon further decrease in gap width, they can be induced to migrate again. Apparently, for intermediate gap widths, the upper wall represents a perturbation inducing a shape change towards the third dimension (i.e., the thickness of the cell) that interferes with the in-plane dimension of the cell and renders it less motile. In turn, sufficiently thin slits squeeze the cell, and they become more motile than on a single flat substrate: they can move for lower values of the actin-associated driving since the strong confinement enlarges the region of adhesion, the lamellipodium and hence force generation. This can be further exemplified by studying cells on a substrate below their threshold of motion.
Slowly approaching a second wall sets them into motion (see Supplementary Fig. 1).
We also studied briefly the possibility of guiding cells by spatially modulated slit widths, as shown in Fig. 1c. There, two different cells, one of lamellipodium-type (green) and one of wallpushing type (red; both weakly driven), were steadily moving and adhering to both plates, before encountering a smooth variation in the distance of the walls, orthogonal to the direction of motion. Interestingly, the lamellipodium-type cell is deflected to thinner gaps while the wall-pushing one prefers wider gaps. The former is easily understandable: the part of the cell in the narrower region broadens the active front while the one in the wider region reduces, hence the cell turns. The latter is more subtle: it appears that the wall-pushing cell is somewhat pinned by the gap gradient and that the broader lamellipodium on the thinner side runs around this pinning point, leading effectively to a turn towards the broader gap width, see Supplementary Movie 4 for details.
Effects of substrate curvature and cylindrical confinement. We next quantified the motion of cells on curved substrates exhibiting a single curvature direction (i.e., cylindrical symmetry), including confinement in a tube and motion on a fiber. We later also briefly discuss effects of Gaussian curvature. The signed substrate curvature c = 1/R (with R the signed radius of curvature) is defined as sketched in Fig. 3a: a cell moving on the outside of a cylinder surface along the main axis corresponds to negative and a cell moving inside a cylindrical tube to positive curvature. Figure 3b displays snapshots of the respective cell shapes, clearly showing that they are severely perturbed by the imposed curvature. This is also reflected in the cell's center-of-mass velocity along the cylinder axis, as quantified in Fig. 3c: Inside a tube (right part), cells move faster than on the flat substrate. In contrast, on the outer cylinder surface they move slower and finally get stopped ("stalled") for a critical cylinder radius that depends on the driving force, i.e., on the acto-myosin parameters.
Overall, cell speed increases with signed substrate curvature, passing by the flat state. The linear dependence of the migration speed V on the substrate curvature 1/R can be anticipated from the following geometric arguments. Let us assume that the migration speed V 0 on the flat substrate depends linearly on the width of the active lamellipodium front W, since the propulsion force depends on the total number of pushing actin filaments. The height h of the cell does not change upon slightly curving the substrate, as verified numerically. Assuming further that the length of the cell (in the direction of motion) remains approximately constant, volume conservation then implies that the cell's cross-sectional area (normal to the direction of motion) does not change. Let us approximate for simplicity the cell's cross-section by a rectangle, S = Wh. On a substrate with positive radius of curvature R, i.e., for a cell inside a channel, this rectangle transforms into a cylindrical shell segment of cross-sectioñ S ¼W=R½R 2 =2 À ðR À hÞ 2 =2 ¼Wh½1 À h=ð2RÞ. Using S ¼S, and solving forW for R ) 1 we obtain W ¼ W½1 þ h=ð2RÞ and hence the migration speed can be given as This relation is in excellent agreement with the simulations, as exemplified in Fig. 3c for the weakly driven lamellipodium-type cell by the black line, using the measured velocity and height of the cell on the flat substrate.
For sufficiently high positive curvature, corresponding to thin capillaries, cells adhere and form lamellipodia all around. This corresponds to the cylindrical confinement inside a microchannel, restricting the motion to quasi-1D as studied e.g. in refs. 15,16 . In the case of lamellipodium-type cells, at a certain positive radius of curvature a competition starts between the perturbed crescentshape states, as observed for lower curvature, and bullet-like shapes with cylindrical symmetry. This competition induces more complex dynamics, as discussed below. Finally, for even thinner tubes, the bullet-shape wins and the cell's speed is increasing again with signed curvature. Wall-pushing cells also switch from crescent-like to bullet shape, close to the transition, however, they typically become arrested.
In the regime of complex dynamics marked in Fig. 3c, the cell apparently cannot decide which shape to adoptthe perturbed crescent or the bullet. Consequently, the cells do not move on average, but as they remain active they induce shape and position oscillations inside the tube. In Fig. 4a the center of mass position in z-direction together with a measure for the axial length of the cell were tracked. The former oscillates, corresponding to alternating crescent shapes attached to the top and the bottom of the channel. In between these extremes, however, the two wings of the crescent cells become unstable and fuse to a more bullet shape, giving rise to double period oscillations of the axial extension of the cell. This can be understood from the shapes during half a period shown in Fig. 4b. The detailed dynamics is sensitive to the specific initial conditions. The fact that the weaker driven wall-pushing cell did not display oscillations hints at the existence of a threshold value of activity for the phenomenon. For longer times even more complex oscillations can emerge, involving corkscrew-type deformations, see Fig. 4c.
Motion on fibers. Coming back to negative signed curvatures, we also studied cell motion on fibers, as experimentally studied for instance in refs. 27,28 . As already shown in Fig. 3, cells can move on fibers with radii larger than their unperturbed size. They can even wrap around sufficiently thin fibers, see Fig. 1f, and move if perturbed by adding a small initial polarization. Interestingly, we found a generic alignment effect: cells initially polarized obliquely to the fiber axis tend to align their direction of motion with the axis. To obtain quantitative insight, we studied the alignment time τ as a function of the fiber radius: cells were first allowed to spread on fibers of different radii and then motion was initiated with a substantial off-axis component (typically 50% of the final axial velocity). Figure 5 displays results for wall pushing-type cells. The alignment rate (inverse alignment time; obtained by exponential fits to the off-axis velocity component) is roughly linear in fiber curvature for the currently accessible range of curvaturesnote that a strict validation down to c → 0 would need excessively large simulation boxes. For the highest curvature, c = 1/R = −0.06, the found value of τ~800 roughly corresponds to the time the cell needs to move 20-30 times its own size.  Effects of Gaussian curvature. We next studied the effect of a finite Gaussian curvature c G ¼ 1 i.e., of substrates displaying two non-zero principal curvatures 1/R 1 , 1/R 2 . The simplest case is a sphere (having R 1 = R 2 and c G > 0), whereon cells are found to circle around with constant velocity. A more interesting case is the "sphere-with-skirt" structure recently proposed in ref. 26 , which we studied by depositing cells on the top of the structure (c G > 0, Fig. 6a), on the skirt (c G < 0, Fig. 6b), as well as on the transition region (Fig. 6c). Experimentally it was found that cells can not spread on the top, an effect attributed to stress fibers, which are stiff fibers constituted of acto-myosin. Instead they were found to spread preferentially in the region where c G changes from positive to negative and to circle around the structure. The first effect is not captured by our model, due to the neglect of modeling the stiff stress fibers, i.e., cells deposited on top spread as shown in Fig. 6a and prefer to stay there upon perturbations. However, cells indeed prefer the region where c G changes and then preferentially move along the azimuthal direction, see Fig. 6c. If deposited on the region with c G < 0, cells typically are repelled and move away from the structure as shown in Fig. 6b. Nevertheless, if sufficiently polarized and heading towards the structure, they can be "captured" by and then circle around it.
To sum up the curvature effects, they can all be understood along the line of arguments that led to Eq. (4). Namely, with the signed substrate curvature the cell-substrate interaction increases and concomitantly also the cell's actin polarization field is enhanced. In other words, cells spread more easily on positive curvatures (the inside of capillaries) than on negatively curved substrates (the outside of fibers), and hence the direction of motion aligns in the direction of the maximum signed curvature.
Specifically, for cells moving on the outside of fibers, the maximum curvature is directed along the principle axis (as the non-zero curvature along the circumference is negative) and hence motile cells align along that axis. For cells inside a capillary, on the other hand, the maximum curvature is along the circumference and we indeed find a preference of cells to circle around the cylinder wall. This may also provide insight into the complex dynamics regime, occurring for intermediate curvatures as shown in Figs. 3c and 4. There, two active fronts are circling In the case of lamellipodium-type cells, at a certain positive curvature, a competition between the perturbed crescent-shape and a bullet-like shape with cylindrical symmetry sets in, leading to complex dynamics, namely shape and position oscillations inside the tube. Finally the bullet shape stabilizes for even thinner tubes. Wall-pushing (and hence less driven) cells become typically stopped in this region of competition, before they attain similar bullet shapes at higher substrate curvatures against each other and if they "collide" head-on, shape oscillations are observed as shown in Fig. 4. Note that the level of driving is important again, because the (comparatively) weaker driven, wall-pushing type cells simply stop, not exhibiting shape oscillations. This is to some extent analogous to the behavior in vertical confinement, see Fig. 2b, where the two motile branches for strongly driven cells exhibit bistability while weaker driven cells are stopped. Finally, in the "sphere-with-skirt" geometry, the Gaussian curvature varies from positive (spherical cap) to negative (skirt). The skirt has one positive curvature, the other principal curvature (with respect to the symmetry axis) is negative and increases to zero with the distance from the spherical cap. Hence, there is an increase of mean curvature in the direction of the widening skirt (pointing away from the structure), explaining why cells generally move away from such structures.
Guiding cells by substrate topography. The 3D modeling framework also allows to study cell motion on essentially flat substrates exhibiting simple topographical features, complementing experimental and theoretical analyses of cell guidance on microprinted adhesive patterns 42,59,60 . Although there are many possible surface features, we focused here on a type that has been extensively studied experimentally [21][22][23][24][25] , namely periodic grooves/ ridges on the substrate's surface. Experimentally, such substrates are supposed to model the sub-cellular topography of extracellular matrix fibers, often serving as directional cues for moving cells 61 .
The surface pattern was created by allowing a square wave in the substrate's phase field to equilibrate for a few time steps before fixing the resulting smoothened pattern. The amplitude and the wavelength of the pattern were varied, ranging from small to similar compared to the height and planar extension of the cell, respectively. Cells were initialized at a 45°angle to the periodic grooves and then characterized by their attained direction of motion as shown in Fig. 7a.
We found that long wavelength patterns (but still of subcellular size) led to perfect alignment of motion of both lamellipodium-and wall pushing-type cells along the grooves, see the red circles in Fig. 7a. In addition, cells elongated on the pattern in the direction of motion, see Fig. 7b. The alignment was found to be rather insensitive to the (sub-cellular) pattern amplitude, and it was typically fast, i.e., the cell only needed to move few times its size. In contrast, patterns of smaller wavelengths led to practically no alignment for both cell types and all studied amplitudes (black squares), and also the shapes were less affected. Interestingly, in between there is a wavelength region (blue triangles) where wall pushing-type cells align, albeit very slowly, perpendicularly to the grooves as shown in Fig. 7c. Concluding, in general, cells prefer to align along grooves comparable to their size. They are "attracted" by the grooves and rather "dislike" the barriers.

Discussion
The framework developed and analyzed here is a generic physical model for lamellipodium-based motion. It lacks, e.g., specific regulation mechanisms as well as the cell nucleus (which can influence motion in confinement 8 ) and adhesion is taken into account only on the level of substrate friction. Moreover, experimentally observed migration behaviors depend on the specific cell type and even the phenotype. Direct comparisons to experiments must hence be taken with caution. It is nevertheless constructive to confront the general trends found in the generic framework to available experimental evidence. Cells confined between two co-planar walls have been investigated in refs. 17,18 . However, the focus was on the transitions from mesenchymal (i.e. adhesive, lamellipodial) to amoeboidal (non-adhesive) motion. The transitions were induced by reducing the adhesiveness in a confined surrounding, favoring the latter motility mode. A detailed study of lamellipodium-based motion in varying confinement at constant adhesiveness is lacking to date. The non-trivial results obtained already in our simple physical model, namely bistability and guidance (see Figs. 1c and 2) strongly suggest that it is worthwhile to carry out such a study.
Cells confined in microchannels have been investigated in refs. 15,16 , albeit with a focus on invasiveness and integrin-myosin interplay, respectively. In both experiments, the channel cross-section was rectangular due to microfabrication restrictions. Nevertheless, cells adhered and formed lamellipodia all around the channel 16 , as in our model. Concerning cell speeds, a velocity increase by 30% in the channel compared to the flat substrate was found for wild type CHO cells 16 . The increase was much more pronounced when α4/paxillin binding was negatively regulated, resulting in amoeboidal motion which is, however, beyond the scope of our model. Amoeboidal motion was probably also induced by the confinement in ref. 15 studying pancreatic cancer cells, where a velocity increase by almost a factor of three was found. The generic model suggests a typical velocity increase by 40-50%, see Fig. 3c. This increase is a consequence of the confinement-induced shape change and the concomitant reorganization of the active lamellipodial front. Consequently, the experimentally observed increase could, at least in part, be traced back to the physical interaction between the environment and the cell and its internal actin structure.
Cell motion on fibers of various sizes has been experimentally studied for instance in refs. [27][28][29] . The early study in ref. 27 investigated various cell types on fibronectin covered glass slides compared to thin fibronectin fibers deposited on glass slides. Little influence on cell speed was found, but the persistence of motion was strongly increased by the presence of the fibers. This is confirmed by very recent experiments 29  fibroblasts (albeit without tracking the motion explicitly) on highly and weakly curved cylinders. There it was found that cells align more efficiently on cylinders with radii of the order of their own size compared to larger radii. These results are in agreement with the alignment effect found in the model and shown in Fig. 5. Collective cell migration on fibers has been also studied 28 .
Although different from single cell motion, an arrest of motion for submicron fibers was observed. Since at this scale the front consists of a single cell, this arrest may be related to the one observed in the model (see Fig. 3c) for small fiber radii. The effect of Gaussian curvature, see Fig. 6, was already critically discussed above in view of the recent experiments of ref. 26 . For a better reproduction of the experimentally observed behavior, the model has to be augmented by actin stress fibers and possibly the cell's nucleus.
Cell motion on periodic grooves or ridges on the surface of planar substrates has been extensively studied experimentally [21][22][23][24][25] . The alignment along the grooves as obtained in our model for patterns with long but still sub-cellular wavelengths (see red circles in Fig. 7a) has been found in all these experiments. In addition, the elongation of cells along the direction of motion (see Fig. 7b) was found in refs. 24,25 for both T cells and Dictyostelium cells. Nevertheless, the perpendicular motion could also be observed implying less elongated cells 25 , as recovered by the model (Fig. 7c). Ref. 25 even determined a "contact guidance efficiency" that first increased with the pattern wavelength before decreasing again when ridges become too distant for a single cell to sense. Although the model displays fast alignment (within a few cell sizes traveled) for long but sub-cellular wavelengths, which is the optimum range obtained in ref. 25 , the time scale is quite sensitive to the actin parameters (especially θ). Hence, on the modeling side, the model should be tuned to a specific cell type to become quantitative. On the experimental side, it would be very interesting to include systematic variations of the pattern's amplitude, to finally obtain a full picture of the contact guidance mechanism and its determining factors.
Let us also briefly compare our modeling framework to previous 3D cell migration modeling efforts. Previously, a two-fluid framework 37 had been developed that allowed to reproduce shapes reminiscent of both fibroblasts and keratocytes 36 . Starting from a disc-shaped cell one had, however, to assume a certain fraction of the cell's circumference to be (and always remain) protrusive. A more recent approach 1 used active gel theory, but also had to assume a predefined leading edge. Hence while both models are more detailed concerning the internal cellular dynamics, the cell symmetry was broken (and hence cell polarity induced) by construction. Contrary, in our model, the equations of motion are completely symmetric and symmetry breaking emerges spontaneously in a self-organized way, as also recently confirmed by the occurrence of shape waves in a 2D version of the model 44 . Recently, a finite-element framework similar to ref. 36 , but including adhesion, was used to model pseudopodbased migration 38 . All these works only investigated motion on flat substrates, but possibly would in principle allow for a detailed analysis of motion in confinement, curved substrates, and topographical features, as performed here, which would be interesting to critically assess our findings. Concerning motion in 3D fiber matrices, there are already a few modeling studies on ameboidal motion. In ref. 62 , a hybrid agent-based/finite-element model was developed that includes pseudopod formation and even blebbing (via variable actin-plasma membrane linkage), while ref. 39 used a phase field approach using multiple internal density fields and especially an actin activator located but diffusing on the cell's surface. These works are clearly complementary to our work, as they focus on activated pseudopods instead of well-defined lamellipodia. It would be very interesting to study within our model cell motion on an ensemble of few fibers, to compare the two distinct motility modes in this biologically relevant setting.
To conclude, we have demonstrated that the developed threedimensional computational model is capable to describe lamellipodium-driven crawling motion of cells inin principlearbitrarily shaped surroundings and for a variety of cell types. Additionally, in contrast to earlier 3D models 1,36-38 , the approach does not need to pre-define protrusive regions in the cell, as it inherently displays bistability between immobile vs. moving states 2 . We investigated in depth several well-defined and experimentally studied scenarios of geometrical perturbation during cellular motion: such as a systematic variation of the substrate's curvature (from cells on thin fibers to the movement inside a capillary), vertical confinement between two plates, as well as motion on topographically structured substrates. The discovered modulations in cell behavior and the resulting guiding principles of motile cells are due to purely physical effects (confinement, curvature) and their interplay with the self-organized shape and acto-myosin machinery. This study will therefore be fruitful in helping to differentiate such simpler effects from more complex cellular response due to biochemical cues and regulation.
The model was intentionally kept simple. We have already shown for the two-dimensional case that the modeling approach is modular 63 , hence explicit adhesion dynamics and substrate deformation 42,43 , global feedback on actin due to membrane tension 48 , as well as regulation (e.g., via Rho GTPases 64,65 ) or additional components involved in the motility 66 can be added without too much complication. Our approach may pave the way to study other 3D modes of motion beyond lamellipodia in the future, which are often triggered by confinement 8,67 . For instance, a confinement-induced transition of slow mesenchymal cells to a fast ameboidal motion has been found in the slit geometry 17 . Modes of 3D cellular motion distinct from lamellipodia also comprise blebbing 68 , osmotic force 69 (i.e., flow of water through the cell), as well as pseudo-/ lobopodial motion 32,39,62 or the formation of invadopodia 70 . The latter is especially interesting as they allow cells to modify their surroundings and hence feedback on the confinement 8,61 not only by deforming, but also by chemically degrading the extracellular matrix 71 . Another route is related to the cell nucleus affecting migration, especially when the confinement or radius of curvature becomes comparable to the nucleus size. The nucleus can be included into the computational framework via an additional passive phase-field, as done in 2D in ref. 72 . In fact, including a deformable nucleus for motion in complex 3D environments would be a worthwhile endeavor for future investigations. Ultimately, motion inside 3D networks should become accessible 73 , similar as developed in ref. 39 for ameboidal motion employing pseudopods. Also for lamellipodial motion, this situation is arguably the most abundant and generic scenario for 3D motility. There, the mesh size vs. the size and deformability of the nucleus present new limiting steps on motion 8 , including feedback between motility and the genetic program of cells.

Methods
Numerical method and model validation. The Eqs. (1) and (2) are solved on a cubic, typically 256 3 , grid with an operator split Fourier pseudo-spectral code on graphical processing units (GPUs) using the CUDA programming language. The static substrate fields, Φ and Ψ are predefined beforehand. The simulation data were rendered and displayed using ParaView, a 3D data analysis and visualization application.
To validate the model, we studied the spreading dynamics, see Supplementary  Fig. 2, and the onset of motion on flat substrates, see Supplementary Fig. 3. The latter reproduces the 2D behavior described in ref. 2 , see also Fig. 1a). More details on the model implementation and testing can be found in Supplementary Note 1.

Data availability
Data that support the findings of this study are available upon reasonable request from the corresponding author.

Code availability
The code used to generate the data shown here is available upon reasonable request from the corresponding author.