Spontaneous polarization and cell guidance on asymmetric nanotopography

Asymmetric nanotopography with sub-cellular dimensions has recently demonstrated the ability to provide a unidirectional bias in cell migration. The details of this guidance depend on the type of cell studied and the design of the nanotopography. This behavior is not yet well understood, so there is a need for a predictive description of cell migration on such nanotopography that captures both the initiation of migration, and the way cell migration evolves. Here, we employ a three-dimensional, physics-based model to study cell guidance on asymmetric nanosawteeth. In agreement with experimental data, our model predicts that asymmetric sawteeth lead to spontaneous motion. Our model demonstrates that the nanosawteeth induce a unidirectional bias in guidance direction that is dependent upon actin polymerization rate and sawtooth dimensions. Motivated by this model, an analysis of previously reported experimental data indicates that the degree of guidance by asymmetric nanosawteeth increases with the cell velocity. Asymmetric nanotopography biases unidirectional cell migration, yet the underpinning molecular determinants are still unclear. The authors use a three-dimensional model to demonstrate that nanosawteeth induce actin-mediated migration directionality, which is dependent on the cell velocity.

M ost cells produce motion through the coupling of the actin cytoskeleton and the cell membrane to the surrounding substrate. In vivo, cell guidance is affected by a variety of external and internal stimuli, including the extracellular matrix (ECM) 1 , surface rigidity 2 , chemical gradients 3 , and nanotopography. Local nanotopography has a substantial effect on cell behavior. For instance, it has been reported that nanotopographic substrates can lead to cell guidance and can influence cell morphology [4][5][6][7] . Nanoridges can lead to the bidirectional contact guidance of cells, causing cells to speed up and elongate parallel to ridges 8 . In vivo environments are often asymmetric 9,10 , in the form of the porous ECM. Such asymmetries are not often considered in modeling cell migration. Asymmetric nanotopographic substrates have recently been shown to lead to a unidirectional bias in cell migration, with a preferred guidance direction that depends on both the details of the nanotopgraphy 11 and cell line 12 .
Physical modeling has been highly successful in describing cell migration, allowing for the creation of conceptually simple, whole-cell models with wide-reaching implications [13][14][15][16][17][18][19][20][21][22][23][24] . Much of the previous modeling of cell migration has focused on cells on flat surfaces. However, the contact-guidance experiments discussed above suggest that the topography encountered in 3D, in vivo environments can have a profound effect on cell migration. We have previously developed a model that successfully captures some aspects of migration on nanotopographic substrates 25 . Here we extend this model to rows of asymmetric nanosawteeth. The cell is modeled as a deformable boundary in which actin polymerizes at the front, forming lamellipodia and pushing the cell boundary outward. This extended model reproduces the key qualitative aspects of cell motility on asymmetric nanotopography.
Our work reveals that on asymmetric, nanoscale sawteeth, cells exhibit spontaneous onset of motion and directed guidance. For a cell in a deterministic model to move on a flat substrate, the cell must initially be polarized 13,16,19,25 . However, when the same cell model is implemented on asymmetric nanotopography, spontaneous cell polarization is observed. In agreement with previous experimental results 11,12 , our model exhibits biased unidirectional motion on asymmetric nanosawteeth, and exhibits a guidance direction that is dependent upon both the details of the nanotopography and actomyosin dynamics. Our model reproduces qualitative cell-shape phenotypes observed experimentally, including cell elongation parallel and perpendicular to the guidance direction, depending upon the parameters of the nanotopography and the cells. The model also predicts that the degree of guidance by asymmetric sawteeth depends on cell velocity. This prediction is supported by a re-analysis of previous experimental data.

Results
Phase-field model. To explore cell migration on nanotopography, we extended a 3D phase-field model that was originally developed to model cells in confined environments 25 . This model uses two dynamical variables to describe the state of the cell. The first, ρ(r, t), is a scalar phase-field that describes the cell boundary and assumes values between 0 and 1. The cell boundary lies at ρ = 1/2, with ρ > 1/2 inside the cell and ρ < 1/2 outside of the cell. The second is a 3D vector field p(r, t). The direction of p gives the orientation of the actin filaments and the magnitude of p gives the ordering of the filaments. As in the previous 3D phase-field model 25 , we describe the substrate using two auxiliary fields: Φ(r) is a steric field that models the volume exclusion of the cell and the substrate and Ψ(r) restricts actin generation to occur close to the surface. The modeling framework is described by two coupled equations for ρ and p, as detailed in the "Methods" section. In the equations, the unit of length is 1 μm and the unit of time is 10 s.
Here we study the effect of three model parameters on guided cell migration and cell shape. The first, β, is the actin polymerization rate, which determines how fast actin polymerizes at the surface of the cell. The second, θ, determines the preferred angle of actin polymerization with respect to the local substrate geometry. When far away from the substrate, actin polymerizes normal to and away from the local cell membrane. Near the substrate, θ determines the component of actin polymerization normal to the substrate. For θ = 0, actin polymerizes tangentially to the local substrate, whereas for θ = 1 the normal and tangential polymerization rates are equal. The third parameter, σ, determines the rate of acto-myosin contraction. A larger value of σ represents more acto-myosin contraction in the cell. Collectively, these parameters can describe a wide variety of cell migration and shape phenotypes 25 .
Asymmetric nanotopography causes spontaneous polarization. We studied cell dynamics on an array of in-registry rows of asymmetric sawteeth, in which both the sawteeth and the spacings between rows were of subcellular dimensions (Fig. 1a). Such sawteeth have been shown experimentally to be able to bias cell guidance unidirectionally. Sun et al. 11 reported that different nanosawtooth heights lead to opposite guidance directions of Dictyostelium discoideum. Furthermore, Chen et al. 12 reported that the different cell lines were guided in opposite directions on nanosawteeth with the same height. For this study, the sawtooth height was varied between 0.75 μm and 2 μm and the sawtooth length was varied between 2.5 μm and 6 μm. We investigated how this nanotopography affects the onset of cell migration and the corresponding migration and shape phenotype evolution by placing a single cell on the substrate, without initial polarization, and observing the subsequent dynamics. The model produces two major cell-shape phenotypes. The crescent phenotype (Fig. 1b) is representative of cells that are more contractile, such as the M4 cell line that has been studied on asymmetric sawteeth 12 . The wedge phenotype (Fig. 1c) is representative of cells on the lower end of the contractility spectrum, such as D. discoideum, which has also been studied on asymmetric sawteeth 11 . Supplementary Movie 1 shows a model cell transitioning from the crescent to the wedge phenotype. Spontaneous polarization is not observed in phase-field models of cells on flat substrates because the resting state is stable. Thus, the cells in this situation cannot migrate without inducing polarization either externally 26 or via initial conditions 13,27 . However, the asymmetry of the substrates studied here allows cells to generate propulsion even in the absence of an initial polarization. The polarization occurs because the cell adheres to the asymmetric surface, which causes a break in the directional symmetry of actin polarization. Figure 1d shows the polarization vector field for a cell that is polarized spontaneously by sawteeth. There are points of strong polarization underneath the cell that has a tangential component up the sawteeth (against the net direction of cell motion). However, there is strong enough polarization at the leading edge of the cell to propel the cell forward. This cell moves spontaneously with θ = 0. Therefore, the impetus for spontaneous motion on sawteeth is actin polymerization tangential to the substrate. In this model, a cell can only move spontaneously down the sawteeth (Fig. 1e) or at an angle to the sawteeth (Fig. 1f).
The direction of cell guidance is influenced by actin polymerization. In experiments with cells on asymmetric nanosawteeth, the guidance direction has been observed to depend on the sawtooth height and length 11 , as well as on the cell type 12 . We next explore why different types of cells can be guided in different directions on the same sawtooth pattern. When inducing an initial polarization up the sawteeth, we observed three possible outcomes in our model: the cell moves up the sawteeth persistently ( Fig. 2a), the cell initially moves up the sawteeth and then turns around to move down the sawteeth (Fig. 2b), or the cell stops. We induce the polarization up the sawteeth to counter the preferred spontaneous polarization down the sawteeth. Trajectories for these three representative behaviors are shown in Fig. 2c.
To probe the selective guidance direction systematically, we created a phase diagram (Fig. 2d) in which we varied β and θ while keeping the sawtooth parameters the same. In all cases, the cell had an initial upward bias. At t = 250, we measured the center-of-mass velocity. This velocity was used to determine the guidance direction shown in Fig. 2d. The data in this figure indicates that as the actin-polymerization rate increases, the cell is able to overcome the nanotopography guidance direction and move down the sawteeth.
We studied the transition in guidance direction by changing β at a fixed θ of 0.1. Figure 2e shows that for a cell with an initial polarization, there is a sharp transition of guidance direction at β = 6.5. Similar transitions occur at other values of θ. The phase transition is not significantly altered if the phase-field simulation mesh size is doubled (See Supplementary Fig. 1). Below β = 6.5, the cell always has the same final (asymptotic) velocity. This guidance velocity is determined by the sawtooth properties, rather than by the model parameters. Additionally, we observe hysteresis in the guidance direction. If the cell is not polarized initially, the cell can only be guided down the sawteeth. This observation indicates that the preferred guidance direction is down the sawteeth. Finally, a cell that reverses direction does not move solely parallel to the sawteeth. We find instead that a cell can turn to move at an oblique angle to the sawtooth orientation.
The substrate parameters can be used to control the cell shape and the guidance direction independently. Our model exhibits a range of stable shape phenotypes that were studied in relation to the guidance direction. To control shape, we modified the contractility of the cell, σ. Contractile forces control the shape of the cell 28 , and are linked closely to cell motility 29 . Here, we investigate the effect that contractility has on the direction of guidance. sawteeth, given an initial polarization in that direction, counter to the preferred direction of spontaneous motion. By examining the dynamics at individual points in the phase diagram (Fig. 3b), we can separate the velocity trajectories into two distinct regions: the first peak and dip in velocity due to the initial polarization and change in shape of the cell, and the ensuing long-term behavior. The initial peak and dip occur for all values of σ, followed by a positive rebound in velocity that occurs at t ≈ 150. The major differences among the trajectories are in the long-time behavior. For σ = 1 there is a small negative acceleration, but the velocity up the sawteeth is persistent. The long-term velocity oscillations can be attributed to two factors: the cell slows down as it crawls up the sawteeth at its leading edge and small cell shape changes cause fluctuations in the center-ofmass velocity. The persistent velocity and the strong contractility of the cell result in elongation perpendicular to the direction of motion and lead to a persistent crescent shape (Fig. 3c). For σ = 0.5, the long-time positive velocity is not persistent. During the period of the simulation the cell moves up the sawteeth, but with a steadily decreasing velocity at long times. The contractility is high, so the cell is again elongated perpendicular to the direction of motion (Fig. 3d). Finally, for σ = 0.1, the long-time velocity is negative. A short period of positive acceleration at t ≈ 175 is observed but does not lead to a second change of direction. The cell in this case has an oval shape that is elongated parallel to the direction of motion (Fig. 3e).
Comparison to experimental cell shapes and motion. Next, we conduct a re-analysis of previous experimental data of D. discoideum migration on asymmetric sawteeth to capture additional details of cell motion. We find that the deterministic phase-field model can reproduce both the observed cell shapes and the preferred guidance direction. Figure 4a, b shows that D. discoideum can be elongated both parallel and perpendicular to the sawteeth, in analogy to the shape phenotypes seen in the model. As shown in Supplementary Movie 2, D. discoideum cells moving along sawteeth tend to elongate in the direction of motion, in accordance with the less contractile cells shown in Fig. 3. Additionally, in a study of different cell lines migrating on sawteeth by Chen et al. 12 , it was observed that M4 cells, a cancerous mutant of MCF10A epithelial cells, exhibit elongation perpendicular to the direction of motion. This behavior agrees with the observation that cancerous cells are generally more contractile than their progenitor cells 29 . Although asymmetric sawteeth create a unidirectional bias of motion, D. discoideum is capable of moving down the sawteeth (Fig. 4c), up the sawteeth (Fig. 4d), as well as at arbitrary angles relative to the sawteeth.
A variety of cell elongation angles with respect to the sawtooth orientation are observed experimentally, as seen in Figure 4e. The distribution of elongation angles has two local maxima at 0 ∘ (parallel to sawteeth) and 90 ∘ (perpendicular to sawteeth). Fig. 4f shows a combined spider plot of cell trajectories on sawteeth of height 1-2.4 μm. By changing the direction of initial cell polarization, a model trajectory can be produced that travels in any direction, and therefore the model can reproduce the arbitrary elongation angles (Supplementary Fig. 1) that are present experimentally.
The simulated cell trajectories in the continuum model are deterministic, and so cannot reproduce the stochasticity of experimental trajectories. Nevertheless, we observe that the  Fig. 3 The model parameters can be used to control shape and guidance independently. a A phase diagram showing the guidance direction as a function of the contractility σ and the actin polymerization rate β (θ = 0.1) with initial polarization up the sawteeth. Increasing the contractility leads to a reversal in the guidance direction. b Cell migration behavior for different values of β and σ. For σ = 1, the cell's initial velocity decays rapidly to a more stable final velocity. For σ = 0.5, the cell's initial velocity decay mirrors that for σ = 1, before decaying further at longer times. For σ = 0.1, the cell reverses direction at t ≈ 100 to move down the sawteeth. 3D visualizations of the stable cell-shape phenotype for (c) σ = 1, (d) σ = 0.5, and (e) σ = 0.1. As the contractility σ decreases, the cell becomes more oval shaped.  direction of stable simulation trajectories corresponds to the preferential guidance direction of cells on asymmetric sawteeth. To determine the effect of sawtooth dimensions on cell guidance direction we compared the model to experiments of cells on sawteeth of three different heights: 1 μm, 1.8 μm, and 2.4 μm. The model predicts that cells with higher actin polymerization rates will move up the higher sawteeth (Fig. 5a). Because cells with higher actin polymerization rate tend to have higher velocities 30,31 , one implication of our model is that fastermoving cells will move up sawteeth with smaller heights and down sawteeth with greater heights. Therefore, the model indirectly predicts that cells with higher velocity will move up the sawteeth for lower sawteeth heights and down the sawteeth for higher sawteeth heights.
To determine the guidance direction of cells in experiments, we measured the average velocity of the D. discoideum cells as a function of direction relative to the sawteeth, with the angle 0 ∘ defined as moving up the sawteeth (Fig. 5b-d). For the 1 μm-high sawteeth, the average velocity peaks at 0 ∘ (Fig. 5b). This observation is consistent with the phase diagram (Fig. 5a), in which cells with a higher actin polymerization rate move up the sawteeth. For the 1.8 μm-high and 2.4 μm-high sawteeth the average velocity peaked at 180 ∘ . This result indicates that the D. discoideum cells do not reach the threshold velocity necessary to move up these sawteeth. Our observations suggest that asymmetric sawteeth can be used to guide cells with high velocities selectively, and that the velocity cutoff can be modified by changing the sawtooth height.
A toy model for cell migration on nanotopography. To strengthen the predictive power of the full phase-field simulations it is useful to capture the main observations in a toy model amenable to analytical extrapolation. The extrapolated positions can be validated in both the phase-field simulations and experimental observations. The toy model incorporates two key observations from the full phase-field model: the sub-critical onset of cell motion on flat surfaces and the unidirectional bias of motion by the sawtooth substrate. This toy model allows us to investigate stability of cell guidance. If there is no topography on the substrate, the equations of motion can be cast in the dimensionless form: Here ϵ ≪ 1 is the relaxation parameter, β is the driving parameter related to actin polymerization, and v is the center-of-mass velocity of the cell. For β > 2, Eq. (1) possesses three stable fixed In the presence of asymmetric nanotopography, the equation of motion assumes the form Here x 0 is the unit vector in x-direction. In the case of asymmetric sawteeth, f(x) can be written in the form Here a is the spatial frequency (spatial period L = 1/a) and c is the magnitude of the sawtooth modulation. Note that the sawteeth provide a positive bias, because 〈f(x)〉 > 0. Despite their simplicity, Eqs. (2, 3) capture many salient features exhibited by the fully 3D, continuum phase-field model as a function of the driving parameter β (Fig. 6). In particular, for pure one-dimensional motion (i.e., in the x-direction), for a large enough value of β the model exhibits two stable solutions, with motion up or down the sawteeth (Fig. 6a). For smaller values of β, the model only shows spontaneous motion down the sawteeth. For even smaller values of β, no motion is possible, all in agreement with the phase-field model.
The two-dimensional version of Eqs. (2, 3) yields a surprising prediction that motion up sawteeth is always unstable (Fig. 6b). An analytical stability analysis (Supplementary Note 1) shows that a cell that is initially polarized in the direction up the sawteeth will, after some time, make a U-turn and begin to move down the sawteeth. However, as β increases the reversal time diverges (Fig. 6c). Thus, the model indicates that the rich migration behavior observed both in the phase-field simulations and in the experiments is a synergistic effect of the sawtooth substrate and the cell's intrinsic threshold for motion. By incorporating only three properties, the toy model leads to three nontrivial predictions: hysteresis between motion up and down the sawteeth, the existence of a critical sawtooth height above which motion up sawteeth is impossible, and the instability of motion up the sawteeth. These predictions indicate that the minimal model ingredients needed to capture velocity-based selective guidance are sub-critical onset of motion on flat substrates and unidirectional bias by the sawtooth substrate.

Discussion
We have demonstrated that a 3D phase-field model enables a rigorous investigation of cell motility on asymmetric nanotopography. Previous work has focused on cell-scale or system-scale gradients. Here we demonstrate that the model can also capture how cells integrate subcellular asymmetry into cell scale guidance. The phase-field model captures the spontaneous onset of motion and the experimentally observed range of shapes and motions of cells that sense asymmetric guidance cues that are an order of magnitude smaller than the scale of the cell.
Cell-migration experiments on asymmetric sawteeth have demonstrated that the guidance direction is affected by the sawtooth dimensions 11 and the specific cell line used 12 . The phase-field model reproduces the key phenomena observed in these experiments, including spontaneous polarization, unidirectional guidance, and elongation phenotypes. By tuning model parameters, we can change the guidance direction and cell shape in accordance with the behavior observed for different cell lines. The correspondence between the model and experiments yields insights into the relative importance of actin polymerization, contractility, or other cell properties in promoting cell-typedependent guidance behavior.
On flat surfaces, the phase-field models require an initial polarization 13,25 or a strongly randomized internal actin density 32 to observe directed cell motion. In contrast, asymmetric sawteeth provide guidance cues that cause the cell to move spontaneously. Consequently, surface asymmetries may be one factor that initiates directional migration in vivo and in natural environments.
Our model shows that the actin polymerization rate is a primary driver of the cell guidance direction on asymmetric sawteeth. There is a transition in guidance direction from down the sawteeth to up the sawteeth at a critical value of the actin polymerization rate. This critical value is dependent on other model parameters as well, such as the sawtooth height. A toy model based on findings from phase-field simulations verifies this relationship and gives analytical evidence that motion up-sawtooth is unstable. The actin polymerization rate is directly related to the cell velocity, so faster cells should tend to move up the sawteeth. In agreement with the predictions of our model, experimental data for cell guidance on asymmetric sawteeth show velocitydependent selective guidance. Additionally, the model predicts that increasing the sawtooth height should increase the velocity threshold for cells to move up the sawteeth. We verified this prediction by showing that faster D. discoideum cells are preferentially guided up asymmetric sawteeth that are 1 μm high. However, on asymmetric sawteeth with heights of 1.8 μm and 2.4 μm, faster cells are still guided down the sawteeth.
In our model, contractility affects both the shape and the guidance direction of a cell. Larger values of contractility are associated with the cell moving up the sawteeth. Thus, more contractile cells are not guided in the more typical direction on asymmetric sawteeth. Studies have linked higher contractility to invasiveness in cancer cells 29 , and this connection may be related to the change in guidance direction that we observe. The simple physical model presented here can explain many features of cell motility on asymmetric nanotopography. Because the phase-field model is highly modular, it can be extended further to obtain a fuller picture of cellular dynamics on asymmetric nanotopography. An even more realistic model could include the modulation of specific or non-specific adhesive effects (e.g., focal adhesions). Additionally, the actin dynamics can be modified to include the effect of regulatory proteins via an excitable network model such as local excitation global inhibition biased excitable network 33 or linked excitable networks 34 .
Although our approach captures many important aspects of subcellular environmental asymmetries on cell migration, in vivo the dynamics can be affected by many other factors. It is known that the cells can remodel, and even degrade, the local topography, especially in the context of some invasive cancers 35,36 . Many cancer cells form invadopodia, i.e., actin-rich protrusions of the plasma membrane that are associated with degradation of the extracellular matrix. Incorporation of such phenomena into a phase-field model can shed light on cancer invasiveness and metastasis and is highly desirable. However, a 3D phase-field model on deformable and degradable substrates is a formidable computational task that we leave to future study.

Methods
Phase-field equations. The following model was used to describe the cell migration in 3D: Here, the phase-field ρ describes the cell shape, p characterizes the actin orientation, Φ(r) models a steric field moderating the cell's interaction with the surface, and Ψ(r) restricts actin generation close to the surface. The entire motility mechanism is captured by Eqs. (4,5). Actin polymerizes near the cell boundary and produces a protrusion force via the polymerization ratchet mechanism 37 . As the cell moves forward, the actin polarization at the back of the cell degrades based on a model of acto-myosin contraction 13 . The cell volume is conserved as a fixed point through the term δ[ρ]. This model of cell motility has been shown to reproduce the lamellipodium-based motion of keratocytes on flat surfaces 13 . We use a modified version of the model that accounts for cell surface tension 38 . Surface tension is introduced using the exponential term e −TcΔS , where T is the surface tension strength, c is the curvature, and ΔS is the difference between surface area at a given time step and a reference surface area. This term makes the actin polymerization rate lower at points of high curvature, thereby preventing the cell from tearing apart. We introduced this term by necessity based on the high curvature of the nanotopography studied here.
Numerical method. The model used here is an extension of a previously published method 25 . Equations (4,5) were solved on a 100 μm × 100 μm × 25 μm rectangular grid (corresponding to 576 × 576 × 128 mesh points) in the domain using a split operator method on a periodic domain. The algorithm was implemented on graphical processing units using the NVIDIA CUDA programming language. The diffusion terms were calculated in Fourier space, and other operators were calculated using finite-difference methods. The curvature c in Eq. (4) was calculated using c ¼ ∇ Á ∇ρ j∇ρj . Unless otherwise mentioned, the model parameters used were those given in Table 1.
The substrate with asymmetric sawteeth was implemented via prescribing the static fields Φ and Ψ. The substrate studied for this paper was a sawtooth of length 4.2 μm and a default height 1.4 μm; other heights were used as discussed above.
Cell tracking. Cell tracking was performed on bright-field movies of D. discoideum moving on asymmetric sawteeth reported in a previous paper 11 . First, the threshold of the fast Fourier transform (FFT) of the images was used to remove the bright spikes associated with the periodic sawteeth. A low-pass Hamming filter was applied to the FFT. An inverse FFT was then applied to recover an image with sawteeth removed. Next, the image contrast was changed so that 1% of the brightest and darkest parts of the image were saturated. A morphological top-hat filter was applied to remove high-frequency noise. Finally, the image was binarized using the Otsu thresholding method 39 . The cells were tracked using the MATLAB regionprops function to find elongation and centroid location. The cell trajectories were found using a particle tracking algorithm fed the centroids of the cell locations.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

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

Code availability
Code that created the data for this study is available upon reasonable request to the corresponding author.