Optogenetic switching of migration of contractile cells

Cell crawling on flat substrates is based on intracellular flows of the actin cytoskeleton that are driven by both actin polymerization at the front and myosin contractility at the back. The new experimental tool of optogenetics makes it possible to spatially control contraction and thereby possibly also cell migration. Here we analyze this situation theoretically using a one-dimensional active gel model in which the excluded volume interactions of myosin and their aggregation into minifilaments is modeled by a supercritical van der Waals fluid. This physically simple and transparent, but nonlinear and thermodynamically rigorous model predicts bistability between sessile and motile solutions. We then show that one can switch between these two states at realistic parameter ranges via optogenetic activation or inhibition of contractility, in agreement with recent experiments. We also predict the required activation strengths and initiation times.

Crawling cell migration is an ubiquitous process in animal tissue and plays a crucial role in e.g.development, wound healing, immune response and cancer metastasis [1].In addition, an increased interest in synthetic cellular systems drives the need for understanding the minimal components required for cell motility [2].The most essential element of cell migration is the symmetry break between front and back: while the front uses actin polymerization to push the membrane forward, the back uses myosin II contractility to pull the rear forward.Traditionally, these processes are considered to be coordinated by gradients in biochemical activity, most prominently the antagonistic signaling pathways of Rac/Cdc42 and RhoA for front and back, respectively [3].
A striking feature of locomoting cells is the bistability of their motility behavior: sessile cells can be stimulated into migration by the application of physical stimuli if the applied stimulus is sufficiently large to polarize their cytoskeleton [4,5].Very recently, it was shown that the direction of cell crawling in channels can be reversed by optogenetic stimulation effectively decreasing myosin II contraction at the back [6].Because this situation is essentially one-dimensional, here the bistability is mediated by polarization, in contrast to cell migration on two-dimensional substrates, where bistability can also result from cell shape changes [5,7].Surprisingly, the optogenetic experiments also revealed that high myosin II activity by itself can prevent reorientation [6].In general, optogenetic perturbations of cell contractility have revealed that contractile cells usually do not work at saturation, but at an intermediate setpoint of tension that allows for up-and downregulation [6,8,9].These observations shed new light on a long-standing question in the fundamental understanding of cell migration, namely how cell migration works and can be controlled in purely contractile cells.
FIG. 1. Cell crawling depends strongly on the spatial distribution of myosin motors.A homogeneous distribution does not induce intracellular flows and therefore corresponds to a sessile cell (top).In contrast, a myosin gradient corresponds to a motile cell (bottom).Optogenetic activation of contractility (red) can be used to switch between these two states.
The natural framework to understand cytoskeletal flow within cells and the role of contraction is active gel theory [10,11], which has been used early on to describe cell migration [12].However, obtaining bistability and switching is difficult within this framework.Previous attempts have relied on introducing nonlinear saturation terms [13,14], but the recent optogenetic experiments [6] demonstrate that the system is not in saturation yet.We also note that the earlier work used an ideal gas description for the myosins, which does not reflect its high density and its property to aggregate into large complexes.
Here we go beyond these assumptions and propose to describe the myosins as a supercritical van der Waals (vdW) fluid, a concept suggested before for other protein systems [15].In the myosin context, it accounts both for the crowdedness of the cytosol and the aggregation of myosin II into so-called minifilaments, which are supra-molecular clusters that lead to persistent contraction of the actin cytoskeleton.The vdW approach results in nonlinearities but is consistent concerning linear irreversible thermodynamics, since the driving force, the gradient in chemical potential, is still in linear order.Here we show that this description provides a physically reasonable and transparent framework to explain the experimentally observed bistability and to predict the effect of optogenetic control (Fig. 1).We also parametrize our model and demonstrate that our predictions both agree with the recent experiments on motility reorientation [6].Finally we use our new theory to predict different experimental protocols to achieve experimental control of cell migration.
To model the cell we assume the constitutive relation of an infinitely compressible active gel with a linear dependence of the active stress on the myosin concentration field c(x, t): Here v(x, t) is the flow velocity field, η the shear viscosity, σ(x, t) the total stress field and χ the contractility per motor protein.We assume viscous drag with the substrate, ∂ x σ = ξv, with a friction coefficient ξ.The cell is considered to have a variable length, with left edge l − (t) and right edge l + (t), and an elastic boundary condition σ = −k(L − L 0 )/L 0 , where L = l + − l − is the cell length and L 0 its reference length.We finally assume that the flow fields at the boundaries agree with cell movement (v(l ± ) = l± ).This implies that we disregard the effect of actin polymerization, which could be incorporated as a modified boundary condition [16].
To determine the dynamic equation for the myosin concentration, let us start with the chemical potential µ c of the vdW fluid [17] Here N A is the Avogadro number, k B the Boltzmann constant, T temperature, b the vdW excluded volume and a the average value of the attractive interaction energy per unit concentration.λ th = (2π 2 /mk B T ) is the thermal wavelength.According to linear irreversible thermodynamics [18], the diffusive particle flux J D follows from the gradient of the chemical potential, J D ∝ ∂ x µ.The continuity equation ∂ t c = −∂ x J D then yields a diffusion equation with a concentration-dependent diffusion coefficient where we identified the saturation concentration c ∞ = 1/N A b and the attractive energy e A = 2aN A /k B T .D(c) has a singularity for c = c ∞ , reminiscent of hard core repulsion, increasing diffusion at high concentrations.
The energetic term leads to a slow down of diffusion for intermediate concentrations due to the attraction.We limit our discussion to the supercritical vdW fluid, i.e., the temperature is above the critical temperature k B T c = 8a/27b [17], corresponding to attractive energies e A < e (c) A = 27/4c ∞ .We note that the effects of this concentration-dependent diffusion coefficient have been observed experimentally [19], for instance in the context of binary liquids [20] and colloidal suspensions of hard spheres [21].It has been used in models for bacterial growth [22] and very recently also for excluded volume effects of myosins in cells [23].
In addition to performing diffusion, the myosin motors are advected by the active gel, ).We consider no-flux boundary conditions at the edges, i.e., ∂ x c(l ± ) = 0. Following earlier work along these lines [13,16], we non-dimensionalize length by L 0 , time by L 2 0 /D, stress by k, and concentration by c 0 = c dx/L 0 .We then map the problem on the interval [0, 1] using u = (x − l − )/L and get three dimensionless model parameters: the hydrodynamic damping length L = η/(ξL 2 0 ) arising from the competition between viscous and frictional dissipation; the Péclet number Pe = k/ξD describing the importance of advection versus diffusion; and myosin contractility P = χc 0 /k.The inverse Péclet number A = 1/Pe can also be interpreted as adhesion strength, because a large value of A corresponds to strong friction if L 2 A = Dη/kL 2 0 is kept fixed, which is a combination of quantities that typically cannot be changed in experiments.In the following, adhesiveness A and contractility P are considered as the main parameters, as experimentally they are known to control transitions in cell state.
Defining the cell center G = (l + + l − )/2 and the advection velocity v(u) = − Ġ + L(1/2 − u), and rescaling c = Lc, we finally arrive at our central equations, namely the following boundary value problem (BVP): with the boundary conditions σ(l ± ) = −(L − 1) and ∂ x c(l ± ) = 0.As we will show below, this system can be comprehensively analyzed using a combination of analytical and numerical methods.
The parameters of the model can be estimated from experimental data for crawling cells.Following earlier work [16], one obtains L 2 = 1.25 and P = 0.1.A can be determined from k = 10 4 Pa [13,16,24,25], ξ = 2 • 10 14 Pa s /m 2 [16,26] and D = 0.7 • 10 −12 m 2 / s to be A ≈ 1/70.A rough estimate for the volume of one (unclustered) myosin motor is 10 2 nm 2 • 100 nm and the myosin concentration in cells is of the order of c 0 ≈ µM [27].This implies the estimate c ∞ = 100, not accounting, however, for crowding in the cell or finite thickness of the cortex, which should decrease this number.Therefore here we use c ∞ = 10, which implies e [17], with e min (in units of k B T ) the binding energy in the Lennard-Jones potential of the vdW fluid.
We will now study the effect of the nonlinear diffusion and show that it results in bistability between sessile and motile solutions.To find the steady states one assumes l± = V , with velocity V , L = 0 and steady profiles, i.e., σ and c only depend on u.One then obtains two coupled ordinary differential equations.The case D(c) ≡ 1 has already been studied in [28,29]: two non-motile solution families exist with flat stress profiles σ ≡ −(L − 1) and lengths L± = (1 ± √ 1 − 4P)/2.More complex solutions bifurcate from these branches.However, only two were found to be asymptotically stable, the trivial branch L+ and a motile branch bifurcating from it, with a peak in the motor concentration at the trailing edge [28].
For nonlinear diffusion, a flat stress profile is still a steady state solution.In the following we focus on the bifurcation from the sessile, flat-stress state (L = L+ , σ = 1 − L+ , ĉ = 1/ L+ ) to the first motile state, as our numerical results suggest that these are again the only stable solutions for experimentally relevant parameters.Fig. 2(a) shows the continuation in Pe of the solutions' cell length and velocity.We see that approaching the critical e (c) A from below renders the supercritical pitchfork bifurcation toward the motile solution to be subcritical, implying bistability.Analytically obtained bifurcation points (cf.SI), indicated as circles, agree with the numerics, showing that increasing the attraction e A in addition decreases the value of Pe at which the bifurcation occurs.Hence attractive interactions both induce bistablity and reduce the motility threshold.We also investigated the full, time-dependent BVP numerically, using the discontinuous Galerkin finite elements method [30,31] (see SI for details), and found that indeed both solutions marked in Fig. 2(a) as solid curves are stable in the bistable regime (cf.SI, Fig. S1).
Using advanced continuation methods (branch point and fold continuation [32]), we determined the boundaries of the three different regimes -sessile, bistable and motile -as shown by the state diagram in Fig. 2(b).Note that we now kept L 2 A fixed as explained above.In Fig. 2(b) we focused on the range around the experimentally reasonable values P ≈ 0.1 and A ≈ 1/77 (a full diagram can be found in the SI).Starting in the bistable regime, increasing adhesion leads to a transition to the non-motile state, while increasing contractility favors the motile regime.Already a change in contractility of 5% allows for these transitions.Note that changing e A shifts all the boundaries and hence, depending on e A , also the opposite scenario can occur, i.e. decreasing adhesion inducing motility.
Finally, Fig. 2(c) shows the normalized motor concentrations Lc(u) of the three possible states in the coexistence region.The unstable branch displays enrichment of motors at the trailing edge, while the motile branch develops a boundary layer.Importantly, the volume exclusion of the vdW limits the height of the concentration (and stress) peak at the edges even for small A (cf. SI Fig. S1(e)).This has to be contrasted to the linear model, where unrealistically large peaks develop [13], and is in agreement with experimental findings of moderate myosin enrichment [33].The flow profiles u shown in Fig. 2(d) indicate the flow to the trailing edge.The flow velocity is maximum in the motor-enriched boundary region.This, together with the attraction and subsequent minimum diffusion, promotes the formation of the myosin layer.
We stress that for e A = 0 (no attraction, only excluded volume, i.e., the Tonks gas) the concentration peaks are still limited, but bistability does not occur.Thus, accelerating diffusion at high concentrations, as commonly assumed [23,28], is not sufficient in this thermodynamically consistent model to achieve bistability.We also find that one needs to be careful using a truncated Taylor expansion as this may result in changes in the criticality of the bifurcation for certain (unphysical) parameters of the full excluded volume interaction (see SI for details).
Having established that the model displays bistability for experimentally realistic parameters, we next ask if optogenetics can be used to switch between the sessile and motile solutions.Optogenetic control of contractility usually exploits light-induced recruitment of a GTP-exchange factor to the cell membrane, which in turn activates the RhoA signaling pathway and thus leads to an increase in myosin II contractility [8,9].An alternative way has been described very recently: there a chemotaxis receptor was optogenetically activated in neutrophils, which promotes Rac/Cdc42-activity and decreases contractility.Both optogenetic strategies can be implemented in our model by introducing an optogenetic contribution to contraction, P → P + EΞ, with a shape function Ξ encoding the spatio-temporal activation and E the activation strength [16].E is positive (negative), depending on whether one activates (inhibits) myosin II contractility.We consider a box-shaped function within an activation region U act , i.e., Ξ(u, t) = 1 only if u ∈ U act and t ∈ [t on , t off ] with turn-on/turn-off times t on and t off .
Fig. 3 shows that both activation in the front half (a) and inhibition in the back half (b) can be used to induce reversals of direction of cell migration.While the first protocol has not been experimentally explored yet, both the simulated trajectories and changes in length for the latter correspond well to the recent experimental results [6].Note that initiating motility through activation would also speak in favor of contractility saturation not being central in this context.In the activation protocol the cell contracts as contractility is increased.The effect of the perturbation builds up throughout the activation period, since length and velocity are governed by the integrated active stress (see SI).For the inhibition scenario the effect is opposite, as we inhibit in the half with higher initial concentration: an immediate length response and a more gradual velocity change is obtained.In particular, we predict a decrease of |V | after turn-off, exactly as observed experimentally [6].
We next address the question if optogenetic control of contractility can be used to initiate or arrest motility of cells inside the bistable regime.In Fig. 4(a) we started with the non-motile steady state and activated the cell's left half.Motility is indeed initiated for E as small as 1.2% of P. Increasing E further leads to faster initiation.For the smaller E shown, the perturbation is not sufficiently strong for the induced flow to overcome diffusion and the system cannot leave the basin of attraction of the sessile state.In Fig. 4(b) we started with the motile steady state (moving to the right) and activated the right (leading) half.Compared to the case in (a), now larger perturbations are required, because in the motile regime the advection from motility is dominating and has to be overcome.Arrest is possible only when fine-tuning the turn-off of the optogenetic signal: it has to occur in the "re-symmetrized region" that belongs to the basin of attraction of the sessile solution; activating beyond this point rather induces reorientation.Again, larger strengths E lead to faster arrest/reorientation. Having demonstrated the possibility to initiate or arrest cell migration by optogenetics, we finally predict the corresponding time scales.In Fig. 5 we show the time t ini at which the steady state velocity is reached.We find that t ini is larger for stronger adhesion, as one would expect.Initiation is faster for larger activation strengths, with an asymptotic dependence of t ini ∼ (E/P) −1 , cf.Fig. 5(a).Using continuation of the optogenetically perturbed system (cf.SI), we determined the lower initiation boundary for different adhesions, see Fig. 5(b).We find that increasing adhesion not only affects the stability of the steady states, but also slows down the dynamics and increases the necessary activation strength.Note that this can be tested experimentally, as activation strength has been shown to depend on the illuminated area in optogenetic experiments [8].
This The bifurcation points of steady state solutions of the nonlinear system [Eq.( 4) in the main manuscript] from the trivial branches can be shown to be given by solutions of with ] and ĉ = 1/ L± .This implicit formula, generalizing a criterion given in [1] to the nonlinear diffusion case, implies that the bifurcation structure stays similar to the linear case D(c) ≡ 1. Figure S2 shows the relaxation dynamics of the model, determining the stability of the branches.The full bifurcation diagram can be obtained using the continuation method [2] and is shown in Fig. S3.

Description of the used Finite Element Method
The full, time dependent BVP [Eq.(4) in the main manuscript] was solved with the discontinuous Galerkin (dG) finite elements method, accounting for the heterogeneous diffusion constant via the symmetric weighted interior penalty scheme for the description of the diffusion term [3] and for advection via upwinding.Time was directly discretized and the time-dependent equations for c integrated with an implicit Euler scheme.Length and cell center position were integrated via explicit Euler stepping.The solver was implemented with FEniCS [4].

Diffusion of Tonks gas and comparison to Taylor expansion
For the Tonks gas of hard spheres on a line [5] the nonlinear diffusion coefficient is equivalent to the van der Waals diffusion coefficient without any attractive energy (e A = 0), i.e., analogously to Eq. (3) in the main text, (S2) * Corresponding author: schwarz@thphys.uni-heidelberg.de where c ∞ = 1/N A b with hard sphere diameter b.This model, henceforth called the Tonks gas, has the second order Taylor approximation Note that the linear term in the van der Waals fluid with attractive energy is (2/c ∞ − e A )c, which means that a diffusion coefficient with vanishing linear term, as suggested in Ref. [6], is equivalent to e A = 2/c ∞ , which is still supercritical and thus falls within our framework.Consequently, neglecting a linear term means that attraction has implicitly been considered.Figure S1 depicts the bifurcation diagrams for the Tonks gas model, its second order Taylor approximation, the quadratic model with e A = 2/c ∞ and consistent volume exclusion, and its Taylor approximation, the latter suggested in Ref. [6], for different saturation concentrations c ∞ .Note that in the Taylor approximations we have no divergence at c = c ∞ and thus can look at saturation concentrations below the steady state concentration c ∞ < ĉ = 1/ L+ ≈ 1.1.For the Taylor approximation of the Tonks model we find bistability only in this unphysical regime, while including thermodynamic consistent volume exclusion does not yield bistability.For the truncated quadratic model bistability is generally also in the unphysical region or very close to c ∞ = 1/ L+ , not consistent with experiments.
Effectively, truncating the Taylor expansion means that we introduce higher order terms cancelling the higher order terms of volume exclusion.These terms are not consistent with volume exclusion and cannot be neglected as in the nonlinear regime large concentrations are observed [7].

Length dynamics governed by integrated active stress
To see that length and velocity dynamics are governed by the integrated active stress and hence the integrated optogenetic signals, the equation for stress [Eq.( 4) in the main manuscript] .Note that we assumed the signal Ξ to only depend on the internal position u, i.e., the cell is in a non-motile steady state or the activation region is moved along with the cell.Similarly to Ref. [8] we find that the steady state length is dominated by the active stress term, i.e., the last integral of the concentration field c in Eq. (S5).
Using the Green's function for σ from Eq. (S4), one can see that the velocity V only depends on the antisymmetric part of the integrated active stress −(P + EΞ)c with respect to the cell center, u = 1/2, weighted by an appropriate integration kernel from the homogeneous solution.Thus the velocity also depends on the integrated active stress.

Continuation with optogenetic activation
To study the effect of optogenetic activation we included optogenetic activation in the active stress in the equation for stress [Eq.( 4 Pe/1000 Full quadratic

AFIG. 2 .
FIG. 2. Bistability of cell migration.Panel (a) gives the cell velocities V and the length differences from the sessile state, ∆L = L − L+, for the obtained solution branches as a function of Péclet number Pe for different (supercritical) attractive energies eA.The bifurcation points are marked with circles.Stable (unstable) solutions are shown as solid (dashed).Panel (b) shows the state diagram for eA = 0.63 and for L 2 A = 1.25/77 fixed (containing parameters, which cannot be changed by the cell, cf.text).Depending on adhesion strength A and contractility P, one finds a sessile, bistable or motile regime.Parameter values estimated from experiments are marked with a circle and used in panel (c).The solid/dotdashed curves correspond to the loci of the pitchfork/saddle node bifurcation.Panel (c) shows normalized motor concentration profiles for experimental parameters (and V ≥ 0) in the bistable regime for the stable motile, the unstable motile and stable sessile solutions, as solid curves.The flow velocities are shown in dashed.

FIG. 3 .
FIG. 3. Reversal of cell migration by optogenetics.(a) Positive activation strength E = 0.07 at the front and (b) negative activation strength E = −0.07 at the back both lead to persistent reversal of migration.The upper panels show kymographs, i.e. material cell points as function of time in lab coordinates; myosin concentration is color-coded.The lower panels show cell velocity V and length L scaled by initial length L(0).Time periods of activation/inhibition are shaded.

FIG. 4 .
FIG. 4. Optogenetic initiation/arrest of motility.(a) Motility can be initiated by optogenetically perturbing the flat motor concentration in (the left) half of the cell; activation E/P = 1.2%, 0.8% for A, B. (b) Motility can be arrested or reoriented when a motile steady state is activated in the leading (right) half; activation E/P = 6.7%, 6.6%, 4% for A, B, and C. Solid curves are simulations and dashed (dotted) lines represent stable (unstable) steady states; motile (non-motile) states are in grey (black).Time periods of activation are shaded.

FIG. 5 .
FIG.5.Initiation times of motility in the bistable regime.(a) For different levels of adhesion the initiation time tini for motility differs, where A = 0.013 corresponds to Pe = 77.Time tini decreases for larger activations E/P, asymptotically decaying as tini ∼ (E/P) −1 .(b) For larger adhesion tini grows larger, with concomitantly increasing minimal activation E/P necessary for motility initiation.The dashed line is the initiation threshold, obtained using continuation.
S4) arXiv:2206.05915v1[q-bio.SC] 13 Jun 2022 is integrated.The steady state length including optogenetic activation is then given by Lact E Ξ(u)] c(u)du, (S5) with the stress deviation field from the boundary condition s(u) = σ(u) + (L − 1) ) in the main manuscript], −P/L c → −(P + EΞ)/L c, with a sharp continuous representation of the left-half box function Ξcontinuation to calculate the bifurcation diagram for different activation strenghts E for the same parameters as in Fig.2(a) in the main manuscript with e A = 0.63, see Fig.S4(a).For optogenetic activation the pitchfork bifurcation separates into two saddle-node bifurcations, where the loss of stability of the previously non-motile solution is marked by the movement of the saddle node-bifuraction beyond the value of Pe.In Fig.S4(b) we continue the non-motile solutions (V = 0) for different Pe in E. The saddle-node bifurcation, where branch switching occurs, if activated beyond, marks the activation threshold for motility initiation for fixed Pe.The stability threshold, shown in Fig.5in the main manuscript is determined as the loci of this bifurcation in the (A, E) plane for fixed L 2 A = 1.25/77.

FIG. S1 .
FIG. S1.Bifurcation diagrams in nondimensionalized velocity V and length L as functions of the Péclet number Pe for different nonlinear diffusion models with volume exclusion.(a) Tonks gas model, which is the van der Waals model without attraction.(b) Second order Taylor approximatione of the Tonks gas.(c) Full quadratic model, i.e., the van der Waals fluid with eA = 2/c∞, such that the linear term cancels.(d) Second order Taylor approximation of the quadratic model.The term linear in concentration c vanishes.In the truncated models we find bistability only in an unphysical regime, see text.Parameters: P = 0.1, L 2 = 1.25 research was conducted within the Max Planck School Matter to Life supported by the German Federal Ministry of Education and Research (BMBF) in collaboration with the Max Planck Society.We also acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Projektnummer 390978043.USS is member of the Interdisciplinary Center for Scientific Computing (IWR) at Heidelberg.