*Article* **Cell-Substrate Patterns Driven by Curvature-Sensitive Actin Polymerization: Waves and Podosomes**

#### **Moshe Naoz and Nir S. Gov \***

Department of Chemical and Biological Physics, Weizmann Institute of Science, P.O.B. 26, Rehovot 76100, Israel; moshe.naoz@gmail.com

**\*** Correspondence: nir.gov@weizmann.ac.il

Received: 16 February 2020; Accepted: 17 March 2020; Published: 23 March 2020

**Abstract:** Cells adhered to an external solid substrate are observed to exhibit rich dynamics of actin structures on the basal membrane, which are distinct from those observed on the dorsal (free) membrane. Here we explore the dynamics of curved membrane proteins, or protein complexes, that recruit actin polymerization when the membrane is confined by the solid substrate. Such curved proteins can induce the spontaneous formation of membrane protrusions on the dorsal side of cells. However, on the basal side of the cells, such protrusions can only extend as far as the solid substrate and this constraint can convert such protrusions into propagating wave-like structures. We also demonstrate that adhesion molecules can stabilize localized protrusions that resemble some features of podosomes. This coupling of curvature and actin forces may underlie the differences in the observed actin-membrane dynamics between the basal and dorsal sides of adhered cells.

**Keywords:** actin waves; curved proteins; dynamic instability; podosomes

#### **1. Introduction**

The actin cortex of cells is the prominent driver of membrane shape deformations, which exhibit a huge variability, from propagating waves to stable protrusions. It is often observed that the actin-membrane dynamics of adhered cells is very different between the basal and dorsal sides. One of the major differences between the two sides is that on the basal side the membrane is held at close proximity to the solid substrate (in the range of 10–100 nm, depending on the stage of the adhesion [1]), while on the dorsal side the membrane is usually free to deform into the surrounding fluid. In this paper we explore theoretically the actin-membrane dynamics in the presence of the confinement of the substrate, when the actin polymerization is nucleated by curved membrane complexes.

Cells exhibit a variety of propagating waves of actin polymerization on their basal plasma membrane, which are observed under many conditions such as the initial formation of adhesion [2] and during cell motility [3–5]. When these waves propagate on the dorsal side of an adhered cell, or along its perimeter edge, they are accompanied by large membrane deformations. However, when these waves propagate along the basal membrane, at the interface between the cell and the underlying solid substrate, such membrane deformations have not been unambiguously observed. These basal actin waves have been studied intensively [6–9], and many of their features exposed. Mostly these waves have been treated in the framework of reaction-diffusion models [9], where membrane deformations do not play a role.

In previous works [10,11] we have investigated theoretically and experimentally the possible role of curved activators of actin polymerization in the propagation of membrane-actin waves. In these works the positive feedback is in the form of an actin nucleator that has a convex shape (such as the I-BAR protein IRsp53 for example [12,13]), such that it tends to accumulate at the tips of membrane protrusions that are driven by the actin polymerization force. The negative feedback, which is necessary for wave propagation, can be provided by the contractile force of myosin-II [10] or the recruitment of concave-shaped actin nucleators (such as the BAR family proteins, Tuba for example [14]) [11]. More recent work proposed that the negative feedback for propagating basal actin waves arises from the actin network itself [15].

In this paper, we explore the dynamics of the membrane-actin system in the presence of only the convex nucleator, but in the presence of a confining boundary which represents the effect of the solid substrate. When there is no confinement, our model predicts that this system can become unstable and drive the spontaneous initiation of membrane protrusions through a Turing-type instability [16–20], as is indeed observed in experiments [21–24]. We show that in the presence of a confining boundary this system indeed supports protrusions, which are however modified compared to those growing on a free membrane: protrusions may split, and may even convert into propagating rings. These theoretical results may explain some puzzling features observed for actin waves that propagate at the substrate-attached cell surface, such as their tendency to form doublets of concentric actin fronts [9,25,26].

In the last section we demonstrate that by adding adhesion of the membrane to the substrate localized protrusive structures can be stabilized, and these share some features with localized adhesion structures called podosomes.

#### **2. Results**

#### *2.1. Expanding Ring of Membrane-Actin Wave*

Our model is based on the description of the membrane shape in terms of a single height variable *h*(*x*, *y*), which is appropriate for small membrane deformations (Monge gauge). This is applicable for the present system, where the membrane is adhered to a solid substrate that confines the extent of its normal deflection. on the membrane we consider a density field *n*(*x*, *y*) of "activator" proteins, which are both curved and can recruit the polymerization of actin (such as I-BAR protein IRsp53, for example). These "proteins" may therefore represent bound complexes that contain several proteins that together have this combination of properties. The curved membrane complexes can diffuse on the fluid membrane, as well as form high density aggregations. The spontaneous curvature of membrane proteins spans a large range, from being almost flat [27], to having radii of curvature of the order of 10 nm [12] (we used a value of 100 nm in this work). Please note that in general the membrane-bound protein complex can have a radius of curvature that is different from that of its individual components.

We solve the equations of motion for the two fields (Equations (10) and (11)) numerically using an explicit finite difference scheme with periodic boundary conditions. We first investigate the response of the system to a single small Gaussian perturbation (see Supplementary Movie S1, Figure 1). We choose values for the parameters of the model such that we are in the unstable regime, and protrusions grow spontaneously. However, the numerical values of these parameters are not fitted to any particular experimental measurement, and are simply chosen to demonstrate the qualitative behavior over realistic length and time scales.

As shown in Figure 1a, the perturbation initially grows into a protrusion with a lateral width of the order of the most unstable wavelength *λ<sup>c</sup>* (Equation (12)). During this growth stage, the protrusive force of the actin locally squeezes the layer of long molecules (glycocalix) that buffers the outer surface of the membrane from the substrate (Figure 8). Due to the positive feedback, the density of activators increases at the tip of the growing protrusion (note that throughout the paper the plots of the "activator density" is with respect to the background, uniform density *n*-*n*0).

**Figure 1.** (**a**–**d**) Numerical integration of Equations (10) and (11) over a period of 3 min, over a membrane segment of size 10 <sup>×</sup> <sup>10</sup> <sup>μ</sup>m2. The parameter values used: *<sup>D</sup>* <sup>=</sup> 0.1 <sup>μ</sup>m2/s, *<sup>H</sup>*¯ <sup>=</sup> <sup>−</sup><sup>10</sup> <sup>μ</sup>m−1, *<sup>n</sup>*<sup>0</sup> <sup>=</sup> <sup>50</sup> <sup>μ</sup>m<sup>−</sup>2, *ns* <sup>=</sup> <sup>300</sup> <sup>μ</sup>m<sup>−</sup>2, *<sup>A</sup>* <sup>=</sup> 3.8 · <sup>10</sup>−<sup>5</sup> kg <sup>μ</sup>m5 sec<sup>−</sup>2, *<sup>κ</sup>* <sup>=</sup> 20 k*B*T, *<sup>σ</sup>* <sup>=</sup> 8.28 · <sup>10</sup>−<sup>5</sup> kg <sup>μ</sup>m4 sec<sup>−</sup>2, <sup>μ</sup> <sup>=</sup> 1.66 · 106 sec <sup>μ</sup>m−<sup>2</sup> kg<sup>−</sup>1, *<sup>h</sup>*wall <sup>=</sup> 0.5 <sup>μ</sup>m. (**a**) Initial growth of the protrusion, prior to contact with the substrate. Please note that throughout the paper the plots of the "activator density" is with respect to the background, uniform density *n*-*n*0. (**b**) the protrusion after it comes into contact with the substrate and the membrane at the tip becomes flat. Consequently, an activator ring is formed at the edge of the membrane protrusion, where there is large positive curvature. (**c**) the membrane at the disk center has retracted back towards the unperturbed position at *h* = 0, and a secondary inner ring of activators begins to form. (**d**) the membrane ring and two activator rings expand further. (**e**) an illustration of the mechanism that drives the expansion of the protrusion radially outwards. When the membrane reaches the flat substrate its curvature diminishes and the activators are then aggregated at the location of the highest curvature—the shoulders. However, since once the activators aggregate they push the membrane against the substrate which results in the flattening of the shoulders and the formation of new shoulders further away from the protrusion center. (**f**) a diagram of the structure of the expanding ring. Marked in green are the regions high in curvature in the radial direction which is similar in magnitude for both the inner and the outer rings. Marked in red is the outer radius curvature in the azimuthal direction which is positive and marked in blue is the inner radius azimuthal curvature which is negative. Also shown is the concentration of activators which aggregate into two rings at the outer and inner radii of the membrane ring. The concentration of the outer activators ring is higher than the concentration at the inner ring due to the different azimuthal curvatures.

Once the protrusion's height exceeds *h*wall, the substrate limits further growth and the membrane shape tends to match the contour of the substrate. If the substrate is sufficiently flat the activators which were aggregated at the tip of the protrusion disperse. Please note that we do not consider at this stage any adhesion to the substrate, so that the activators remain mobile on the membrane even when it is in contact with the substrate. The result is a rolling instability where the activators continually aggregate at shoulders of the protrusion (Figure 1b,e), which are the location of the highest convex curvature. The aggregation of activators increases the protrusive force exerted on the shoulders which are therefore results in the membrane deformation moving radially outwards. The rolling instability results in the protrusion developing into an expanding circular structure.

We emphasize that this model includes only normal deformations of the membrane, so the actin force does not directly push the membrane sideways along the substrate. The movement of the protrusion laterally is driven by the flow of the curved activators, and the coupling to the protrusive force of the actin polymerization. Please note that a similar behavior is expected for curved activators that adsorb in a curvature-dependent manner from the cytoplasm [11].

The membrane's initial shape is an expanding cylinder and the activators form a ring at the membrane perimeter of the protrusion (Figure 1b). Once the radius of the membrane cylinder is sufficiently large, the membrane at the center of the cylinder, which is no longer supported by a surplus of actin protrusive force, falls back to the initial height (at *h* = 0). This happens due to the inherent repulsion between the membrane and the substrate, cause by the "cushion" layer of long molecules (glycocalix) that cover the outer surface of the membrane (Figure 8). When the membrane cylinder changes into a ring shape, a secondary inner ring of activators forms at the inner shoulder of the membrane ring (Figure 1c,d,f), where there is high convex curvature.

The amplitude of the activators density at the inner ring is initially considerably smaller then the outer ring amplitude. The reason for the amplitude difference is the difference in the mean curvature between the outer an inner rings. While the radial curvatures (the curvature along the radial coordinate centered at the protrusion center) are very similar, the azimuthal curvature (the curvature along the angular coordinate), which is of the order of 1/*R* (*R* is the radius of the ring) is positive at the outer radius and negative at the inner radius (Figure 1f). Therefore, at small radii, the difference in the mean curvatures is significant. The higher convex curvature at the outer ring means that the convex activators aggregate there more and the resulting protrusive force exerted on the membrane is stronger. This imbalance results in the continued outwards expansion of the ring. Please note that the inner actin ring does *not* move inwards, since the curved actin activators flow towards increasing mean convex curvature, which decreases if the ring would shift to a smaller radius. Therefore the inner ring is also propagating outwards, at a velocity which is very similar to that of the outer ring, maintaining a roughly constant distance between them.

However, as the ring radius grows larger, the difference between the azimuthal curvatures at the inner and outer rings diminishes and the difference in the amplitudes of the activators rings (and therefore the protrusive force) decreases, which reduces the speed of the ring expansion (Figure 2a,b). In Figure 2c we plot the activators density at the inner and outer rings as function of the local mean curvature, and in Figure 2d we plot the ring velocity as a function of the difference between the density of activators at the inner and outer rings. The plot shows that the velocity is proportional to that difference, i.e., *V*ring ∝ Δ*n* = *n*outer − *n*inner. The results indicate that the ring velocity is indeed proportional to the imbalance in the pushing force of the two actin rings, and explains why it decreases as: *V*ring ∼ 1/*R* (Figure 2b).

**Figure 2.** (**a**) the radius of the outer membrane disk (and later ring) as a function of time (all panels correspond to the simulation shown in Figure 1). (**b**) a log plot of the expansion speed of the outer radius vs the radius. We see that the graph is curved at small radii (where the ring is actually a disk) but as the radius grows (and the ring forms) the graph approaches a straight line indicating the power law relation: *V*ring ∼ 1/*R*. (**c**) the peak values of the differential concentration *n*-*n*<sup>0</sup> at the inner (red) and outer (blue) rings as a function of the local membrane curvature. We see the the concentrations are linear in the curvature, as given by Equation (3). (**d**) the ring outwards velocity as a function of the difference in concentrations between the inner and outer activators ring.

We can quantify these observations by the following calculation: If we hold the membrane shape constant, we can solve the steady state concentration profile of the curved activators that corresponds to that shape. By taking *n*˙ = 0 in Equation (11) and integrating once we get

$$D\nabla n + \frac{\Lambda \mathbf{x} \vec{H}^2}{n\_s^2} n \nabla n - \frac{\Lambda \mathbf{x} \vec{H}}{n\_s} n \nabla^3 h = 0 \tag{1}$$

we then divide by *n* and integrate again and we are left with

$$D\ln\left(n/n\_0\right) + \frac{\Lambda\kappa\dot{H}^2}{n\_s^2}(n - n\_0) - \frac{\Lambda\kappa\dot{H}}{n\_s}\nabla^2 h = 0\tag{2}$$

For large concentrations we can neglect the first term on the left hand side and get

$$n - n\_0 = \frac{H}{n\_s} \nabla^2 h \tag{3}$$

We therefore find that when the dynamics of the activators is faster than the expansion rate of the membrane deformation, so that the activators' concentration is in a quasi-steady state, we get that the activators amplitude is approximately proportional to the local membrane curvature. In Figure 2c we plot the concentration of the inner and outer rings vs. the mean curvature at these locations. The plot shows the concentrations are a linear function of the mean curvature which confirms that for the parameters used in the calculation, the dynamics of the activators are indeed faster than the membrane dynamics, and the result of Equation (3) is valid.

Since in our model the concentration of actin activators is strongly affected by the local membrane curvature, the dynamics of the ring is sensitive to the topography of the substrate. We illustrate this by simulating the dynamics on a substrate that is roughened with random bumps with an average amplitude of a few tens of nanometers (Figure 3a). Using the same set of equations and the same initial conditions as shown in Figure 1a–d, we now get the result shown in Figure 3b,c. The overall qualitative behavior is similar to the case with a smooth surface, i.e., the formation of an expanding ring-like membrane structure with inner and outer rings of activators. However, the shape of the expanding structure is strongly affected by the substrate roughness and did not retain the circular symmetry it started with. The membrane ring also shows short "finger-like" protrusions extruding radially from the perimeter, which are accompanied by very strong aggregation of activators. The inner activators ring does not extend into these deformations. Due to the surface roughness, and the consequent fluctuations in the membrane curvature, the distribution of the activators inside the membrane ring-like structure becomes very inhomogeneous and fragmented.

**Figure 3.** (**a**) Substrate with random roughness, of Gaussian amplitude, with an average amplitude of a few tens of nanometers. (**b**,**c**) Numerical integration over a period of 3 min of a ring expanding over a membrane segment of size 10 <sup>×</sup> <sup>10</sup> <sup>μ</sup>m2. We used the same parameter values as in Figure 1.

Another illustration of the effects of substrate topography is shown in Figure 4 (Supplementary Movie S2), where a single elongated cylindrical ridge protrudes from the otherwise flat surface (along the *x*-axis). As can be seen, when the membrane ring first reaches the tip of the ridge it is curved backwards with respect to its expansion direction (top middle panel). As the membrane wraps around the protruding ridge, this membrane part develops negative curvature and the ring of actin activators breaks up at that point (top right panel). Only when the inner activators ring forms, this tip of the bump becomes favorable and concentrates activators that begin to push the membrane ring backwards towards the point of initiation (bottom-right panel). on the two side of the elevated ridge the original membrane ring propagates faster than on the flat substrate. This is due to the higher concentration of actin

activators, which is caused by the high positive curvature in the sharp corners that the elevated ridge makes with the flat substrate (bottom middle and right panels).

**Figure 4.** Numerical integration of a ring expanding over a membrane that contains a cylindrical ridge along the *x*-axis. From top-left to bottom-right, each panel shows the system at increasing time, with the membrane shape and activator density in the top and bottom parts, respectively. We used the same parameter values as in Figure 1.

The simulations in Figures 3 and 4, demonstrate that the distribution of actin activators (and therefore actin filaments) can become highly fragmented due to surface undulations, while the enveloping membrane structure remains continuous and smooth. It is therefore not straightforward to relate the actin signal to the membrane topography when interpreting experimental images of such cell-substrate waves. It also shows that topographic features can cause local direction reversals and break-up of these actin waves.

#### *2.2. Array of Localized Protrusions*

We now study the conditions that may stabilize localized structures driven by the same curved activator we used so far. In Figure 5 we show that starting with different initial conditions may lead to a quasi-stable array of localized structures. A random noise in the initial membrane height or activators density, instead of a single perturbation, gives rise to multiple protrusions, each with a lateral size comparable to *λ<sup>c</sup>* (Figure 5). When these protrusions reach the substrate they undergo the same process of flattening and expanding that we saw before for a single protrusion. During their expansion, the distance between these protrusions naturally decreases. The interaction between these protrusions is repulsive due to the positive membrane curvature region trapped between neighboring protrusions, similar to those observed in other curved membrane-bound aggregates [28]. Therefore, once they have expanded and reached a distance of order *λ<sup>c</sup>* from each other their expansion is halted, and they stabilize. Over the course of the stabilization, the protrusions expand into all the space that is available by the repulsive interactions between neighbors, leading to the elongation of some of the protrusions. Furthermore, if two protrusions were initially formed in close proximity, the energy barrier between them is decreased and they can coalesce into a single elongated protrusion (arrows in Figure 5b–d point to such a process).

Please note that in these localized protrusions the actin forms small rings within each protrusion, due to the same process we found in the expanding ring (Figure 1), on a smaller scale. These protrusions, which may seem stable, are highly sensitive to small perturbations, since they are stabilized only by their mutual repulsion due to the local membrane-driven barrier. Over long times we expect noise to cause them to shift and coalesce. Non-linear terms drive coalescence over such barriers, over a long time [28].

**Figure 5.** Simulation over a period of 3 min over a membrane segment of size 5 <sup>×</sup> <sup>5</sup> <sup>μ</sup>m2, with Gaussian noise in the initial membrane height with variance of 10nm. The membrane's random initial deformations develop into protrusions of lateral size ∼ *λc*. (**a**) the initial growth period before most of them make contact with the substrate. (**b**,**c**) the stabilization period of the protrusions, which elongate into the available space and may coalesce with very proximal protrusions. (**d**) the final steady state of the system. We used the same parameter values as in Figure 1.

#### *2.3. Adhesion-Stabilized Localized Protrusion: The Podosome*

To stabilize an isolated protrusion on the basal side, using the curvature-actin mechanism that we propose, we need to stabilize the localized actin core and prevent the tendency of the protrusion to expand outwards as a ring (Figure 1). An example of a localized adhesion structure on the basal size of many cell types is the podosome [29]. Podosomes are actin-rich protrusive adhesion structures formed on the membrane of several cell types, and have been implicated in the processes of cell migration, tissue invasion and extracellular matrix (ECM) degradation. The mechanisms that give rise to podosome formation, and their large-scale organization in the cell, are still poorly understood and are the subject of ongoing current research [30]. Podosomes are typically formed in monocytic cells such as macrophages, osteoclasts and dendritic cells and similar structures called invadopodia have been observed in carcinoma cells [31–35]. They are relatively dynamic, formed and destroyed in the span of a few minutes and are formed only on the interface between the cell and a substrate. Furthermore, the I-BAR protein IRsp53 was found in podosomes [36,37], suggesting that the basic ingredients of our model may apply for this cellular structure.

The podosome's actin core is surrounded by an adhesion ring, which we did not include so far in the theoretical model, and we therefore suggest that this component may stabilize the core and prevent its ring-like expansion. We propose that the adhesion molecules form a diffusion barrier that greatly inhibits the diffusion of the membrane-bound actin nucleators (Equation (14)), and thereby stabilizes the localized core. In addition, there is evidence that *concave* curved membrane proteins, such as BAR domain proteins, accumulate at the podosome [38], and can further act as a diffusion barrier [39]. we incorporate the adhesion into the model with the same approach that we used to incorporate the actin, namely all the proteins involved in the adhesion process, from plaque proteins that form a scaffold around the actin core to the integrins which connect between the cellular membrane and external ligands, are grouped into a single component which we denote the "adhesion proteins". We consider adhesion proteins to be membrane proteins that have two possible states: a non-adhered, freely diffusing state with concentration *gf* and an adhered, immobile state with concentration *gb*. The transition from non-adhered to adhered state is only possible when the distance between the membrane and the substrate is small enough, for the membrane-bound integrins to bind to the substrate. Please note that we do not explicitly describe the inter-podosome actin network [34], but rather focus on modeling a single, isolated podosome.

In addition to the geometric constraint, the binding rate of the membrane-bound adhesion proteins depends on the application of tensile forces, as integrins are known to exhibit catch-bond properties [40,41]. In the vicinity of the actin core of the podosome, the tensile force is thought to arise at the outer edge of the core, where the actin filaments flow towards the actin core and apply a pulling and shearing force that facilitates integrin adhesion [42–44] (Figure 6a).

We combine the two properties that affect the adhesion listed above, in a very simplified way, in the following equations for the binding/unbinding rates of the adhesion proteins (used in the first order kinetics Equation (15))

$$k\_{on}^{\mathcal{S}} = \left. k\_{on,0}^{\mathcal{S}} f(h) \right| \nabla n| \tag{4}$$

$$k\_{off}^{\mathcal{S}} \;=\; k\_{off,0}^{\mathcal{S}} \tag{5}$$

where *f*(*h*) has the profile shown in Figure 6b so that binding is only permitted close to the substrate. The last term in Equation (4) describes in the simplest way the fact that the adhesion is dependent on the spatial gradients in the actin force that is applied on the membrane. In addition, we note that the adhesion strength that determines *k g on*,0/*k g off* ,0 is affected by the substrate stiffness and chemical composition.

As before (Figure 1), we investigate the system's behavior due to a small Gaussian perturbation in the membrane height (Figure 7). Numerical integration of Equation (15) shows that the perturbation grows into a protrusion and the curved activators aggregate at the tip of the protrusion. When the protrusion reaches the substrate, the activators disperse from the center towards the shoulders. However, at the same stage, adhesion proteins bind to the substrate around the activators and inhibit their dispersion. If the adhesion proteins aggregate quickly enough they can trap the activators in the protrusion core, despite the negligible curvature that the membrane at the center possess due to the confinement by the flat substrate. However, the activators in the core may still form a small ring due to small changes in membrane curvature (Figure 7). When we choose a shorter binding protein, we get a core of smaller radius, which is uniform, i.e., the activators do not form a ring

shape. We therefore demonstrate that the adhesion ring can indeed function as a diffusion barrier that stabilizes the actin core. On long time scales, where the membrane-bound actin nucleators may detach from the membrane, the actin core may decay, and this process could limit the podosome lifetime.

**Figure 6.** (**a**) Schematic illustration of the actin cables around the core of the podosome, where the treadmilling flow induces shearing and pulling forces near the membrane. These forces are thought to activate integrin-based cell adhesion [30]. (**b**) the function *f*(*h*) (Equation (4)) that we used to limit the binding of the adhesion proteins only to within a distance of 2*l*<sup>0</sup> from the substrate (here *h*wall = 0.5 μm).

Overall, the stable podosome-like structure that our model produces exhibits many properties of the podosome. This serves to demonstrate that a model with very few components can give rise to spontaneous formation of membrane protrusions at the basal side of cells, which form an adhesion complex that closely resembles podosomes. It is sometimes observed that the actin core of podosomes may be slightly depleted at its center near the membrane [45], which is a feature that can also appear in our model (Figure 7).

Please note that our model contains a single type of curved actin "activator", while podosomes seem to possess a complex composition of actin filaments [46], as well as a complex and dynamic inter-podosome actin-myosin network [47]. Future modeling of these cellular structures, could involve more of these components, allowing for the simulation of their large-scale dynamics, where podosomes form macro structures comprised of many podosomes [44,48]. For such long-term simulations we will also need to include processes that limit the lifetime of individual podosomes. In most cells, podosomes are organized in a cluster of uniform distribution with a characteristic distance between, undergoing processes of fusion and fission [49]. This organization and dynamics qualitatively resembles the dynamics shown in Section 2.2 (Figure 5).

**Figure 7.** Numerical integration of Equation (15) for the case of constitutive active activators over a period of 2 min over a membrane segment of size 5×<sup>5</sup> <sup>μ</sup>m2. The parameter values used were the same as used previously, with the following changes: *<sup>A</sup>* <sup>=</sup> 1.9 · <sup>10</sup>−<sup>5</sup> kg·μm5· sec−2, *Dg* <sup>=</sup> 0.1 <sup>μ</sup>2/sec, *k g on*,0 <sup>=</sup> 0.05 <sup>μ</sup>m5/sec, *<sup>k</sup> g off* ,0 <sup>=</sup> <sup>1</sup> <sup>μ</sup>m2/sec, *<sup>l</sup>*<sup>0</sup> <sup>=</sup> 0.04 <sup>μ</sup>m and *<sup>g</sup>*<sup>0</sup> <sup>=</sup> <sup>50</sup> <sup>μ</sup>m−<sup>2</sup> (the initial uniform concentration of the free adhesion proteins). The perturbation develops into a protrusion and when it reaches the wall the activators at the core start expanding into a ring but this expansion is inhibited by the adhesion ring that forms around it. The actin core is then trapped by the adhesion ring and the structure stabilizes.

In several cell types, such as osteoclasts, large collections of podosomes exhibit a transition to an expanding multi-podosome ring structure. The ring is densely populated with podosomes, has a width of a few podosomes and expands at a speed of ∼1–2 μm/min [50,51]. The podosomes in the multi-podosome ring are immobile and the ring's outward expansion is achieved by a treadmilling manner: the podosomes decay at the inner part of the multi-podosome ring and form preferentially at its outer edge [49]. These multi-podosome rings move outward and merge with each other until they reach the cell periphery, where they may stabilize as a podosome belt ("sealing zone") [52]. In certain cells the rings seem to originate at the same locations at roughly regular intervals [51].

Within the current model we cannot account for the detailed dynamics of the podosomes within the expanding ring, but we can speculate that its outwards expansion may be related to the mechanism that drives the expansion of the actin ring described in Section 2.1: When a podosome in the ring decays, its constituent proteins diffuse away, and due to the ring-link deformation of the membrane the curved proteins tend to aggregate more strongly at the outer edge of the ring compared to its inner edge (see Figure 1). This mechanism will therefore naturally increase the rate of podosome formation on the outer edge of the ring, compared to its inner edge, and over time cause the outwards expansion of the podosome ring. This expansion will depend on the rate of podosome initiation and decay that enables this effective "podosome treadmill"y.

#### **3. Discussion**

The results of our model can give a very natural explanation to several puzzling features observed in experiments. One such feature is that actin waves at the cell-substrate interface are observed to expand as a single ring of actin polymerization, but beyond a certain size an inner ring of actin that follows the outer one appears. The inner ring is often weaker than the outer ring, as our model predicts. This feature was first noted by Vicker [25] in *Dictyostelium discoideum* amoebae, and more recently this feature was studied in great detail [9,26].

The actin fronts observed in these experiments are very often broken and fragmented, with numerous breaks appearing along the ring [25,26]. This feature is a natural consequence in our model of the sensitivity of the actin concentration to the surface topography (Figures 3 and 4), due to the curvature sensitivity of the actin nucleators. Our model further predicts that the actin polymerization will tend to concentrate where the substrate has concave corners, as along the sides of elevated ridges (Figure 4). This prediction is in agreement with observations of crawling *Dicty* on surfaces with patterned ridges [53].

To conclude, we have shown some aspects of the dynamics of the actin-membrane system, when driven by convex actin activators and confined by the substrate. We show that the confinement itself provides a source of negative feedback that can drive propagating fronts. These simulations are simplified and contain several assumptions, which may not apply to all the biological cases. Furthermore, we do not claim that the actin waves along the basal membranes of cells are driven solely by the mechanism that we describe. Clearly complex reaction-diffusion feedbacks play an important role in the propagation of actin waves in cells [4,9,25,26,54–56]. Our work may motivate further studies of models that include both the reaction-diffusion dynamics and the membrane shape, coupled by curved membrane complexes that nucleate actin polymerization [15]. Since reaction-diffusion models lead to rich dynamics, as does the curvature-actin coupling [19], we expect that models with both features could open up new classes of cell membrane dynamics to explore.

Regarding localized structures at the basal membrane, our model predicts that these may be stabilized by the formation of adhesion around the actin core, as observed in podosomes. Furthermore, if the membrane can be maintained in a curved shape at the tip of the protrusion (rather than flatten), we predict that the curved activators will be less strongly dispersed (if at all), and the lifetime of the localized protrusion extended. This prediction is in agreement with observations that podosomes preferentially form along grooves where the membrane naturally has the curvature of the protrusion tip [57,58], and on rough surfaces [59]. When the protrusion is able to penetrate into the substrate, as occurs in invadopodia [60], the curvature at the protrusion tip is also maintained and this can also stabilize it. Future studies could explore the membrane dynamics predicted by this model on curved substrates, during phagocytosis for example.

#### **4. Materials and Methods**

#### *4.1. Model Without Membrane-Substrate Adhesion*

The model has two variables, *h*(*r*, *t*), *n*(*r*, *t*), that describe the local height deformation of the membrane from its uniform state, and the local density of the membrane-bound, and curved, activators of actin polymerization, respectively. We will work in the limit of small membrane deformations, which allow us to treat the elastic energy of the membrane due to tension and bending in the quadratic limit [16,17]. This is an approximation which may be justified due to the presence of the confining boundary that naturally limited the amplitude of the membrane deformations. Non-linear effects arise in our calculations only from the conservation of the activator field *n*. For simplicity we do not consider the process of binding and unbinding of the actin activators from the membrane, which can be added in the future [11].

We start with the free energy of the membrane and actin activators [16,17]

$$F(h, n) \quad = \int d^2r \left[ \frac{1}{2} \kappa \left( \nabla^2 h - H \frac{n}{n\_s} \right)^2 + \frac{1}{2} \sigma \left( \nabla h \right)^2 + k\_B T n \ln \left( n \right) \right] \tag{6}$$

where *σ* is the effective membrane tension, *κ* the bending modulus, *H*¯ the spontaneous curvature of the actin activators and *ns* their saturation density. In this treatment the two-dimensional Laplacian ∇2*<sup>h</sup>* is the local mean membrane curvature, and we keep the entropic term only at the lowest order, valid for low protein densities *n ns*.

We model the external barrier (substrate) as a one sided harmonic potential (Figure 8) that affects the membrane if its height coordinate *h* exceeds the barrier height coordinate *h*wall. We also subject the membrane to a similar force (though smaller in magnitude) if its height is lower than the initial height (at *h* = 0) to account for the overall average rigidity of the cortical cytoskeleton. The wall and cytoskeleton interactions can be inserted as potentials to the free energy of the system, in the form

$$V\_{\text{wall}} = \begin{cases} \frac{1}{2} F\_{\text{wall}} \left( h - h\_{\text{wall}} \right)^2 & h > h\_{\text{wall}} \\ 0 & h \le h\_{\text{wall}} \end{cases} \tag{7}$$

$$V\_{\text{cyt}} = \begin{cases} \frac{1}{2} F\_{\text{cyt}} h^2 & h < 0 \\ 0 & h \ge 0 \end{cases} \tag{8}$$

where *F*wall,*F*cyt determine the stiffness of these potentials.

**Figure 8.** Illustration of the model. The substrate ("wall") acts as a spring with a force *F*wall(*h* − *h*wall) (right pink arrow) when the membrane height *h* exceeds the boundary location *h*wall. The cytoskeleton acts in the same way only in the opposite direction (left pink arrow) and is typically softer, i.e., *F*cyt < *F*wall.

The equation of motion for the membrane height is given by [16,61]

$$\dot{h} = -\int \Gamma(r - r') \frac{\delta F}{\delta h} d^2 r \tag{9}$$

where Γ(*r* − *r* ) is the Oseen tensor for the hydrodynamic interactions through the surrounding fluid. We replace this long-range interaction kernel with a local on, as is often used for membranes that are highly confined by the cytoskeleton (and here also by the substrate) [16,17]. We therefore take: Γ(*r* − *r* ) = μ*δ*(*r* − *r* ), where μ is the drag coefficient of the membrane. We can now use Equations (6) and (9) to write the equation of motion for the membrane height

$$h = \mu \left[ -\kappa \nabla^4 h + \frac{\kappa \hat{H}}{n\_\epsilon} \nabla^2 n + \sigma \nabla^2 h + A(n - n\_0) - F\_{\text{wall}} (h - h\_{\text{wall}}) \theta \left[ h - h\_{\text{wall}} \right] - (1 - \theta[h]) F\_{\text{Cyl}} h \right] \tag{10}$$

where we used the step function *θ*[*h*] to implement the truncated forces of the confining potentials of Equations (7) and (8). The fourth term on the r.h.s. denotes the force due to actin polymerization, which is proportional to the density of actin activators with proportionality factor *A*. We subtract the average density of the actin cortex in this term to denote the fact that the cell membrane tends to be pushed away from the substrate by a layer of extracellular molecules (called the glycocalix) [1,62]. These molecules act as molecular cushions and maintain weak (non-specific) cell-substrate adhesion, and maintain osmotic pressure on the cell membrane.

The dynamics of the actin activators is given by [16,17]

$$\dot{m} = \Lambda \nabla \left( n \nabla \left( \frac{\delta F}{\delta h} \right) \right) = D \nabla^2 n + \frac{\Lambda \kappa \ddot{H}^2}{n\_s^2} \nabla \left( n \nabla n \right) - \frac{\Lambda \kappa \ddot{H}}{n\_s} \nabla \left( n \nabla^3 h \right) \tag{11}$$

where Λ is the mobility of the activators in the membrane, and *D* = Λ*kBT* is their diffusion coefficient.

We note that the model presented here considers activators that are permanently bound to the membrane and respond to the curvature by flowing in the membrane. Alternatively, the activators can be considered to adsorb to the membrane from the cytoplasm in a curvature-dependent manner [11].

Linear stability analysis of these equations of motion (Equations (10) and (11)) indicate that the system is unstable to small perturbations in either the membrane shape or activators density, for a range of parameters. For negligible membrane tension the most unstable wavelength is

$$
\lambda\_{\rm c} \simeq 2\pi \sqrt{\frac{k\_B T \mu}{A \dot{H} \mathfrak{m}\_s}} \tag{12}
$$

We chose the parameters to have this wavelength of order 1 μm.

The numerical simulations were done using an explicit finite difference scheme centered in space and forward in time, with periodic boundary conditions. The Cartesian grid used in the simulations was 0.025 × 0.025 μm in size.

#### *4.2. Model with Membrane-Substrate Adhesion*

When including adhesion proteins, the free energy of the system becomes

$$\begin{split} F(n,h) &= \int d\mathbf{x} dy \left[ \frac{1}{2} \kappa \left( \nabla^2 h - \frac{\bar{H}}{n\_s} \mathbf{n} \right)^2 + \frac{1}{2} \sigma (\nabla h)^2 + \mathcal{g}\_b \theta (h - h\_\mathcal{g}) \frac{1}{2} k\_d (h - h\_{\text{wall}} - l\_0)^2 + \mathcal{g}\_b \theta (h - h\_\mathcal{g})^2 \right] \\ &\times k\_B T \left( n \ln(n) + \mathcal{g}\_f \ln(\mathcal{g}\_f) \right) \Bigg| \tag{13}$$

where *l*<sup>0</sup> is the length of the external part of the adhesion protein and *hg* is the minimal membrane height required for adhesion to take place. We take this value to be *hg h*wall − 2*l*<sup>0</sup> to allow for variations in the proteins length or membrane fluctuations.

In addition, we assume that the mobility of the actin activators to decrease in regions that have high concentration of adhered adhesion proteins. The reasoning behind this is that the scaffold of plaque proteins surrounds the actin core and restricts the movement of actin filaments, and forms a "diffusion barrier". If actin activators are attached to actin filaments then they too will be restricted. We therefore take the mobility to decrease with the local density of bound adhesion proteins, as follows

$$
\Lambda = \begin{cases}
\Lambda\_0 - \frac{\mathcal{G}\_b}{\mathcal{G}\_{\max}}, & \text{for } \mathcal{g}\_b \le \mathcal{g}\_{\max} \\
0, & \text{for } \mathcal{g}\_b > \mathcal{g}\_{\max}
\end{cases}
\tag{14}
$$

*Cells* **2020**, *9*, 782

Under these conditions we can write the following equations of motion, including the dynamics of the free and bound adhesion proteins (*gf* ,*gb* respectively)

$$\dot{\eta}h = -\mu \left( \kappa \nabla^4 h + \frac{\kappa H}{n\_\delta} \nabla^2 n + \mu \sigma \nabla^2 h + A(n - n\_0) - \theta (h - h\_\delta) g\_h k\_d (h - h\_{\text{wall}} - l\_0) - \tag{15a}$$

$$\left(\theta(h - h\_{\text{wall}})F\_{\text{wall}}(h - h\_{\text{wall}}) - (1 - \theta(h))F\_{\text{tgt}}h\right) \tag{15b}$$

$$\dot{m} = D\nabla^2 n + \frac{\Lambda \kappa H^2}{n\_s^2} \nabla(n \nabla n) - \frac{\Lambda \kappa H}{n\_s} \nabla(n \nabla^3 h) \tag{15c}$$

$$\mathbf{g}\_b = \theta (\mathbf{h} - \mathbf{h}\_\mathcal{\boldsymbol{g}}) k\_{\rm on}^\mathcal{\mathcal{G}} \mathbf{g}\_f - k\_{\rm off}^\mathcal{\mathcal{G}} \mathbf{g}\_b \tag{15d}$$

$$\mathbf{g}\_f = D\_\mathcal{g} \nabla^2 n\_f - \theta (\mathbf{h} - \mathbf{h}\_\mathcal{g}) k\_{on}^\mathcal{g} \mathbf{g}\_f + k\_{off}^\mathcal{g} \mathbf{g}\_b \tag{15e}$$

where *Dg* is the diffusion coefficient of the free (un-adhered) adhesion proteins in the membrane, and the binding/unbinding rates of the adhesion proteins is given in Equations (4) and (5).

**Supplementary Materials:** The following are available at http://www.mdpi.com/2073-4409/9/3/782/s1, Movie S1: Movie of numerical integration of Equations (10) and (11) showing a single growing protrusion, that upon collision with the flat substrate is converted into an expanding ring. This movie corresponds to Figure 1 of the paper. Movie S2: Movie of numerical integration of Equations (10) and (11) showing a single growing protrusion, that upon collision with the flat substrate is converted into an expanding ring. The ring then encounters a cylindrical protrusion along the x-axis of the substrate. This movie corresponds to Figure 4 of the paper.

**Author Contributions:** Conceptualization, M.N. and N.S.G.; numerical simulations, M.N.; writing, M.N. and N.S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Israel Science Foundation grant number 337/05.

**Acknowledgments:** N.S.G is the incumbent of the Lee and William Abramowitz Professorial Chair of Biophysics. This work is made possible through the historic generosity of the Perlman family. N.S.G. wishes to thank Alessandra Cambi and Carsten Beta for useful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
