Modeling Conﬁned Cell Migration Mediated by Cytoskeleton Dynamics

: Cell migration is an important biological process that has generated increasing interest during the last several years. This process is based on three phases: protrusion at the front end of the cell, de-adhesion at the rear end and contraction of the cell body, all of them coordinated due to the polymerization/depolymerization of certain cytoskeletal proteins. The aim of this work is to present a mathematical model to simulate the actin polymerization/depolymerization process that regulates the final outcome of cell migration process, considering all the above phases, in a particular case: when the cell is confined in a microfluidic channel. Under these specific conditions, cell migration can be approximated by using one-dimensional simulations. We will propose a system of reaction–diffusion equations to simulate the behavior of the cytoskeletal proteins responsible for protrusion and contraction in the cell, coupled with the mechanical response of the cell, computing its deformations and stresses. Furthermore, a numerical procedure is presented in order to simulate the whole process in a moving and deformable domain corresponding to the cell body.


Introduction
Actin polymerization in cells is an essential process in a wide range of biological phenomena.The coordinate polymerization of these proteins controls the correct movement of the cell, regulating its rigidity, activating its migration or taking part in the mitosis process.In cell migration, the actin-myosin interaction and the actin polymerization itself generates mechanisms that favor the cell advance: a migrating cell creates a protrusion in the migration direction by actin polymerization at the leading edge; then, the cell adheres the leading edge to the surface and de-adheres at the rear; finally, the whole cell contracts due to forces generated by the actin-myosin network and it moves forward.A detailed description can be found in [1,2].Cell migration is crucial in the development of biological processes such as chemotaxis, cancer metastasis or tissue regeneration.Understanding how the mechanical factors are able to regulate cell movement offers the possibility of developing new therapies to, for example, control the advance of invasive tumor cells.In vitro works allow for considering cell migration as a process regulated mainly by the actin filaments' polymerization inside the cell (see [3]).As mentioned before, cell migration can be simplified assuming a mechanical cycle of three main phases: protrusion due to actin polymerization, adhesion and contraction due to the actin-myosin interaction.Increasing the knowledge of the actin polymerization process in combination with the contractile capacity of the actin-myosin interactions is crucial to understanding the role of mechanical conditions on the cell migration process (see [4,5]).
Cell migration is the aim of a significant number of works, such as in silico and in vitro perspectives.Several mathematical models were proposed in order to simulate this process, but each of them focused on different mechanisms such as protrusion, adhesion or contraction.Complete reviews of these different models can be found in [6] and more recently in [7].For example, Refs.[4,5,8] considered one-dimensional mathematical models to compute the deformation and stress undergone by the cell, considering a viscoelastic behavior, together with the actin and myosin concentrations.George et al. [9] also presented a two-dimensional model to describe cell deformation and movement, coupling its behavior with biochemical processes.Moreover, the relationship between the cell velocity and traction forces due to the involved phenomena (protrusion, contraction and adhesion) were studied in [10][11][12] for the two-dimensional case.Another approach presented in [13,14] considered the cell as a two-phase continuum material, where one of the phases is the cytosol and the other the viscous, actin network capable of contracting and swelling.Finally, in [15], a three-dimensional finite element model to simulate the cell migration on a surface was studied, focusing on the mechanical aspect of the process.
Moreover, regarding in vitro experiments, during the last several years, researchers have focused their interest on three-dimensional cell migration and on the usage of microfluidic devices to determine the mechanical properties of cells.In these types of devices, cells are confined to the walls, and the migration mechanisms are different from the two-dimensional case.For example, Wilson et al. [2] studied leading edge protrusion during chemotactic migration by microfluidic channels.There, they focused on the actin structures driving protrusion in three-dimensional environments, identifying two distinct actin networks, and they proposed the mechanical interaction between them as begin the responsible for forward protrusion.Moreover, Mak et al. [16] presented a microfluidic device to study the cell migration in response to physical spatial gradients.Regarding the models adapted to these types of devices, Hawkins et al. [17] presented a mechanism of cell motility coupling actin polymerization to geometric confinement.Later, in [18], they completed the previous one-dimensional model including the effect of the contractile activity of the myosin motors.Finally, Aubry et al. [19] presented a two-dimensional mechanical model to simulate the cell migration through a micro-channel in order to describe the mechanical behavior of both the cytoplasm and the nucleus.The importance of studying cell migration in confined micro-environments was stated by Stroka et al. [20], gathering different assays developed to provide information about the mechanisms of cell migration.
As a first approach and in order to obtain a simple model to predict the direction and magnitude of cell migration, we focus our interest on the interstitial migration, that is, the cell movement through microfluidic channels.Due to the complexity of the whole cell migration process when considering a three-dimensional scenario, we propose simplifying the process solving a suitable one-dimensional model.This approach is appropriate when considering interstitial migration in microfluidic devices due to the significant differences between the spatial dimensions in the channels.Therefore, the aim of this work is to obtain a one-dimensional mathematical model to simulate the behavior of the main components regulating cytoskeleton (CSK) dynamics in order to predict the cell migration, the actin and myosin spatio-temporal distribution and the deformations and stresses undergone by the cells.
With this aim in mind, and based on the previous works of Gracheva and Othmer [8] and Larripa and Mogilner [4], we will propose a system of coupled reaction-diffusion equations to model the actin polymerization process together with the myosin contractile effect on the cell.Moreover, actin and myosin concentrations are a key factor in the cell migration process since they directly affect the mechanical deformability/contractility of the cell.In addition, we will couple the previous equations with the mechanical response of the cell, computing its deformation and stress.The cell body is viewed as an elastic material, including the contraction and protrusion effects due to myosin and actin concentrations.The main differences of the model proposed here from those mentioned in the bibliography [4,5,8] are considering and computing not only the bound actin and bound myosin concentrations but also the free actin and free myosin densities.Reaction-diffusion equations including effects such as polymerization, depolymerization and positive feedbacks are used to model the dynamical behavior of the actin and myosin densities.Moreover, the proposed model includes also the cell growth domain during the process and the adhesion and de-adhesion phenomenon with the substrate.In order to carry out the numerical simulations of the proposed model, we will use a suitable spatial-temporal integration scheme that allows us to compute the effect of the actin and myosin densities on the cytoskeleton deformation.
The model presented here is a first and simple approach of the interstitial migration process.The final aim is to implement a simple procedure to model the behavior of a cell in a microfluidic channel, mainly in the bifurcations, in order to obtain a quick tool to predict the cell behavior in these types of devices in the laboratory.
This work is organized as follows.In Section 2, we will describe the mathematical model to simulate the migration process.Moreover, Section 2.2 is devoted to the proposed numerical procedure to compute the deformations and stresses of a cell during interstitial migration.In Section 3, we will present the numerical results obtained in the simulation of the considered cell migration process.Finally, Section 3 collects some conclusions and comments about future work.

Materials and Methods
In this section, we present the mathematical model and the different mechanisms considered to simulate the cell migration process.Moreover, we briefly describe the numerical method used to solve the proposed model.

Mathematical Model
We assume that the cell movement is regulated by the CSK dynamics, the actin and myosin concentrations and the deformations and stresses undergone by the cell.
Figure 1 shows a scheme of the different interactions considered in the model: the evolution of the actin and myosin concentrations is governed by means of reaction-diffusion equations; the cellular deformations and stresses are related by means of an elastic constitutive law, where the cell stiffness depends on the actin concentration; the force contraction generated by the cell is due to the myosin bound to the actin; the cell growth during the process is due to the stress generated at the ends of the cell body; and the adhesion/failure is also regulated by the existent stress at the ends, increasing/decreasing the friction with the substrate.Let Ω 0 = (0, L 0 ) be the initial one-dimensional domain corresponding to the cell geometry with L 0 the cell length (see Figure 2).We assume that the cell occupies the region (r(t), f (t)) at each time instant t, where r(t) and f (t) denote the rear and front ends, respectively; thus, r(0) = 0, f (0) = L 0 and Ω(t) = (r(t), f (t)).The problem we aim to solve is to compute the displacement u = u(x, t) and the stress σ = σ(x, t), suffered by the cell body Ω(t) during the time interval of interest (0, T).Note that u and σ are scalar variables since we focus on the one-dimensional problem.In order to compute these unknowns, the actin (a b , a f ) and myosin (m b , m f ) concentrations are also considered, distinguishing between the active and inactive (bound and free) behaviors.

Mechanical Equilibrium
The quasi-static equilibrium of the cell is governed by the mechanical equilibrium equation where g denotes the internal volume forces acting on the cell body.We assume that these volume forces are only due to the adhesion of the cell with the substrate, that is, g = g(x, t) = β v(x, t), and where v(x, t) denotes the velocity of the cell body and β is the friction coefficient, which measures the strength of adhesion between the cell and the substrate.The larger this coefficient is, the larger the friction that is produced.

Constitutive Law
We consider that the cell is an elastic material; therefore, following [8], we can relate the stress with the resulting deformation through an elastic passive stress and a contractile active stress: The contractile stress is due to the myosin proteins bound to the actin filaments in the cytoskeleton of the cell; therefore, we assume that where m b denotes the concentration of bound myosin proteins and η m is a characteristic myosin stress parameter.The elastic passive stress is described through a proportionality relation between this stress and the cell strain.The elastic modulus E(x, t) is considered depending on the position and time since the mechanical properties of the cell change also with position and time, with a dependence on the actin density.Due to this, the elastic modulus E(x, t) depends on the number of bound actin monomers in the cytoskeleton, that is, where E 0 is a constant value.
Summing up, the mechanical behavior law remains: Gathering Equations ( 1) and (2), the proposed differential equation to solve the displacement of a cell is the following:

Reaction-Diffusion Equations
The myosin behavior is described by two variables: the active and inactive myosin (m b , m f ).The active myosin concentration corresponds to the density of myosin proteins bound to the actin filaments and the inactive myosin concentration is due to the free myosin proteins.Analogously, the actin behavior is divided into the active and inactive concentrations (a b , a f ).The active concentration corresponds to the bound actin monomers conforming the filaments found in the cytoskeleton, and the inactive concentration is the free actin monomers flowing in the cell.In order to compute their concentrations, we propose the following reaction-diffusion equations that describe the behavior of the different densities including the diffusion and advection terms: Here, D ab , D a f , D mb and D m f are diffusion coefficients and F a and F m are reaction functions that define the behavior of the specific concentrations.The advection terms in the equations model the transport of the different species in the cell.In this work, we consider the following reaction functions: The reaction functions F a and F m are chosen in order to reproduce the behavior of actin monomers and myosin molecules in the cell cytoskeleton (see, for example, [21]).In these reaction terms, the bound actin a b plays the role of an active protein, whereas the free actin a f acts as an inactive protein.Moreover, with these forces, free actin a f is an activator of bound actin a b , bound myosin m b is an inhibitor of bound actin a b and free myosin m f is an activator of bound myosin m b .In Equation ( 8), k a and δ are the rates of activation and inactivation of actin, respectively, γ is the rate of maximal feedback strength and K is its saturation parameter, whereas, in Equation ( 9) k mp and k mdp are the activation rate and the depolymerization rate of myosin molecules, respectively.
Finally, we need to complete our mathematical model with adequate boundary conditions.We assume conditions of null flux for actin and myosin concentrations (either free or bound), that is:

External Forces
In order to model the cell movement through microfluidic channels, we consider that the cell body is free of forces in the rear end (r(t)) and a surface force is temporarily applied to the cell membrane in the front end ( f (t)) during a determined time interval [t a , t d ].We assume that this load is similar to that produced during an aspiration experiment.The expression for this boundary condition is the following: where h is given by: being h a (t) a linear pressure function such that h a (t a ) = σ a and h a (t d ) = σ d , with σ a and σ d the activation and deactivation pressure values and t a and t d the activation and deactivation time steps.

Cell Domain Growth
When external stress is applied, a growth in the cell length could happen.The model presented here assumes that this growth occurs only if the cell is mechanically deformed up to a suitable stress threshold, denoted by σ g , that is: In the second case, the cell membrane detaches from the actin cortex, creating a new space where the actin filaments in the cytoskeleton can expand and, therefore, the cell domain grows (see Figure 3).This new space is determined by using the one-dimensional relationship between strain and stress in an elastic material.To describe the behavior of the different species in the cell when this new space is generated (that is, σ( f (t), t) > σ g ), we propose simplifying the equations of the problem, neglecting the convection terms and distinguishing between global variables in the former cell domain, denoted with the superscript g (e.g., a g b ) and local variables in the new domain, denoted with the superscript l (e.g., a l b ).Since bound actin and bound myosin have a low rate of diffusion, we can consider bound actin and local bound myosin as purely local variables, that is, we assume that the evolution of these proteins in the new domain is independent of the former cell domain.Moreover, since free actin and free myosin diffuse quickly, we can consider that their local concentration is adjusted immediately to the global value and free actin and free myosin are considered purely global variables.A similar technique was presented in the work of Walther et al. [22] in order to study the stability of a partial differential equations (PDE) problem when a local pulse is perturbing a uniform steady state.Actually, this condition is similar to what occurs in the border of a cell when a membrane perturbation occurs, leading to a change in the cell body.
Therefore, we propose the following simplified equations: where Fa is an adapted reaction function to take into account that bound myosin is almost null in the new cell domain: Finally, it is reasonable to consider only a narrow growth in the domain that has almost null influence on the former cell domain; thus, we can consider that the total quantity of actin (free and bound)-respectively, myosin (free and bound)-is constant, and we can express a f in terms of a g b respectively, m f in terms of m g b .Thus, we obtain the following four ordinary differential equations: Therefore, when an external force is applied, we solve these equations instead of Equations ( 4)-(7).

Adhesion Process
In order to simulate the adhesion process, if the passive stress due to the elastic response reaches a compression stress threshold σ c , the friction coefficient with the substrate, β, is widely reduced to allow the de-adhesion movement, that is:

Numerical Method
In order to solve Equations ( 3)-( 14), we use a suitable finite difference scheme to discretize and numerically compute the solution of the considered problem.
Given initial conditions and choosing uniform space and time discretization, we solve the system at each time step following four steps: • Firstly, we use the stable Crank-Nicholson scheme to compute the displacement u, and consequently the stress tensor σ at each time step through the corresponding constitutive law.• Secondly, we update the actin and myosin concentrations (bound and free) using a similar scheme.

•
If an external stress is applied, the cell is deformed mechanically.If the growth stress threshold σ g is reached, the membrane is detached from the actin creating new space and Equations ( 16) and ( 17) are solved.• Finally, after updating the variables, we check if the compression stress threshold σ c is reached in order to update the friction coefficient with the substrate β.
Figure 4 shows a sketch of the numerical method.At each time step, we compute the cell displacement and stress together with the actin and myosin concentrations and we check if some threshold is reached in order to update the corresponding parameters.

Results
In this section, we present some numerical results for a cell migration process by solving the model proposed above.Moreover, we carry out a sensitivity analysis in order to justify its applicability.

Influence of Actin and Myosin Densities on Cell Mechanics
In order to evaluate the model, numerical simulations are conducted.The aim of these simulations is to determine the influence of the actin and myosin concentrations on the cell movement and deformation.The values of the parameters used in the model are gathered in Tables 1 and 2. We consider the following initial conditions for displacements and stresses together with actin and myosin concentrations: These initial conditions are chosen to start the simulation from an equilibrium state for the different densities at the beginning of the process.
The results obtained for 10 s of simulation are shown in Figures 5 and 6.
Figure 5 shows the time evolution of the actin and myosin concentrations.The bound actin concentration grows with time until the application of the external force h a (t).Once this force is applied, the actin density accumulates in the front end of the cell, due to the generation of new space to polymerize, and after it starts to decrease to recover a more uniform concentration.Contrarily, the free actin concentration decreases with time until the application of the external force.In addition, the bound and free myosin shows an opposite behavior.The bound myosin grows with time far from the front end and the free myosin decreases with time during the whole simulation.Finally, Figure 6 shows the displacement together with the stress undergone by the cell body during the simulation.Note that the deformation of the cell body is due to the myosin contraction effect during the interval (0, t a ); after applying the surface load h(t), an increase in the deformation is suffered and, consequently, a growth in the cell domain is produced; finally, during the interval (t d , T), the contraction effect is the responsible for less increase in the cell body length.The final length of the cell body is 10.711 µm, which generates an increase of 0.711 µm in 10 s of simulation.
Moreover, Figure 7 shows the evolution with time of the rear and front ends of the cell body.Note that the migration process continues after the deactivation of the external forces but with a slower speed.

Sensitivity Analysis
In order to justify the applicability of the model presented here, it is important to analyze the influence of the main phenomena involved in the cell migration process.Due to this, a sensitivity analysis has been carried out, varying some basic parameters in order to better understand the process.We consider the numerical simulation presented in the previous section as the reference case and we analyze the response of the method when varying the following parameters: • the initial cell length L 0 , increasing its value to 20 µm, • the characteristic myosin stress η m , decreasing its value to 0.1 pN/µm 2 , in order to study the influence of the cell contractility in its movement, • the time interval of application of the surface load, deactivating the load at t d = 2.5 s.
Figure 8 shows the results obtained in these three scenarios, comparing them with the results obtained with the reference case.When increasing the initial cell length, the obtained deformation is analogous to the reference case (the cell contracts up to 2 s, then it grows when the surface load is applied, and finally it suffers from a lower increase due to the contraction effect).When the myosin stress η m is reduced, the cell barely contracts and the growth in the cell length is greater (1.274 µm in 10 s of simulation).Finally, when delaying the deactivation time from t d = 2.1 s to t d = 2.5 s, we obtain similar results to the reference case, but the contraction effect is lower and a growth of 4.658 µm is produced.

Discussion and Conclusions
In this work, we introduce a one-dimensional mathematical model to simulate the cell deformation and movement during interstitial migration, that is, when the cell is confined in a microfluidic channel.This model allows for studying the behavior of the cell and the proteins involved in the process: free and bound actin and free and bound myosin.It takes into account not only the mechanical deformation undergone by the cell body, but also the evolution of the protein concentrations during the migration process.
Nevertheless, some hypothesis and simplifications are needed in order to reproduce the cell migration process.As a first approach, we consider that this process can be simulated by a one-dimensional mathematical model.Moreover, this model does not take into account the interaction between the cell cytoskeleton and the surrounding fluid in the microfluidic channel (see, for example, [23]).In this work, we focus on the mechanical response of the cell body with respect to the walls and we neglect the effect of the surrounding fluid during the whole process.Finally, the cell body is considered as a unique entity, without distinguishing between the nucleus and the cytosol.Some authors [19,24,25] simulate the cell nucleus separately during the migration process in microfluidics, highlighting the key role of the nucleus deformability in order to generate adhesions and, therefore, to migrate.More recently, some studies show that cell migration also happens in confined environments with limited adhesion, not only with low adhesiveness (see [26,27]), but also in the absence of adhesion complexes (see [28]).Nevertheless, the effect of adhesion and friction with the walls of the device is considered a surface load in the model presented here.
Regarding the numerical results obtained with this model, as expected, we observe that the larger concentration of bound actin monomers is placed in the front end, where traction forces are applied.Contrarily, the bound myosin is mainly located at the rear extreme of the body cell.Regarding the mechanical behavior of the cell, we notice that the maximum stress is reached when the surface force is applied, and, next, the cell attains a more relaxed state.Taking into account the experimental results obtained for interstitial migration in microfluidic devices in [2,29], the obtained results prove that the proposed model reproduce qualitatively the expected behavior of a cell in interstitial migration.
The next step in this model will be the inclusion of the fluid-structure interaction between the cell and the surrounding fluid in the microfluidic devices in order to simulate the process of directional choice when a cell migrates in these types of confined channels.

Figure 1 .
Figure 1.Scheme of the mathematical model considering the actin and myosin concentrations, the cellular deformations and stresses, the force contraction generated by the cell and the cell growth during the process.For cytoskeleton (CSK) dynamics, free actin a f is considered an activator of bound actin a b , bound myosin m b an inhibitor of bound actin a b and free myosin m f an activator of bound myosin m b .

Figure 2 .
Figure 2. Scheme of the considered experiment.

Figure 3 .
Figure 3. Scheme of the cell domain growth.

Figure 4 .
Figure 4. Flowchart of the numerical method.

Figure 5 .
Figure5.Actin and myosin concentrations during 10 s of the simulation: (Left) bound (solid) and free (dashed) actin concentrations at the initial time and at time instants t = 0, 2, 5, 10 s; (Right) bound (solid) and free (dashed) myosin concentrations at the initial time and at time instants t = 0, 2, 5, 10 s.

Figure 6 .
Figure 6.Stresses undergone by the cell body during 10 s of simulation.The deformation of the cell body with time is also shown.

Figure 7 .
Figure 7. Evolution of the rear and front ends of the cell body during 10 s of simulation.The cell migration process continues after the deactivation of the external forces.

Figure 8 .
Figure 8. Sensitivity analysis of the proposed model.Time evolution of the rear and front ends of the cell are shown.The solid line corresponds to the reference case introduced in Section 3. (a) reference case and three new scenarios; (b) reference case vs. increased initial cell length; (c) reference case vs. decreased myosin stress; (d) reference case vs. delayed deactivation time.

Table 1 .
Numerical values for the mechanical data and phenomena parameters.

Table 2 .
Numerical values for the parameters of the simulation.