The Influence of Nucleus Mechanics in Modelling Adhesion-independent Cell Migration in Structured and Confined Environments

Recent biological experiments (Lämmermann et al. in Nature 453(7191):51–55, 2008; Reversat et al. in Nature 7813:582–585, 2020; Balzer et al. in ASEB J Off Publ Fed Am Soc Exp Biol 26(10):4045–4056, 2012) have shown that certain types of cells are able to move in structured and confined environments even without the activation of focal adhesion. Focusing on this particular phenomenon and based on previous works (Jankowiak et al. in Math Models Methods Appl Sci 30(03):513–537, 2020), we derive a novel two-dimensional mechanical model, which relies on the following physical ingredients: the asymmetrical renewal of the actin cortex supporting the membrane, resulting in a backward flow of material; the mechanical description of the nuclear membrane and the inner nuclear material; the microtubule network guiding nucleus location; the contact interactions between the cell and the external environment. The resulting fourth order system of partial differential equations is then solved numerically to conduct a study of the qualitative effects of the model parameters, mainly those governing the mechanical properties of the nucleus and the geometry of the confining structure. Coherently with biological observations, we find that cells characterized by a stiff nucleus are unable to migrate in channels that can be crossed by cells with a softer nucleus. Regarding the geometry, cell velocity and ability to migrate are influenced by the width of the channel and the wavelength of the external structure. Even though still preliminary, these results may be potentially useful in determining the physical limit of cell migration in confined environments and in designing scaffolds for tissue engineering. Supplementary Information The online version contains supplementary material available at 10.1007/s11538-023-01187-8.


Introduction
Cell migration on two dimensional (2D) substrates and inside three dimensional (3D) environments plays an essential role in many physiological and pathological processes, including embryonic development, wound healing, immune response, cancer progression and metastasis formation (Wolf and Friedl 2006;Bergert et al. 2015;Trepat et al. 2012).The unconfined motion of cells on 2D extracellular matrix (ECM) is a well-studied process and it is conventionally described by continuous and highly coordinated cyclic processes: the elongation of protrusions at the leading edge driven by actin polymerization, the formation of integrin-mediated focal adhesions (FAs), myosin-mediated contraction and the detachment of the trailing edge (Trepat et al. 2012;Balzer et al. 2012;Abercrombie et al. 1970).This classical description requires that specific transmembrane adhesion proteins (integrins, among others) carry intracellular forces from the cytoskeleton to the substrate to propel the cell forward (Bergert et al. 2015;Rafelski and Theriot 2004;Vicente-Manzanares et al. 2009) and it is therefore called adhesion-dependent migration or integrin-mediated migration.
While this mechanism of motion is well understood, the physical challenges that cells have to face when moving in 3D environments are only now receiving more attention, and recent researches indicate that in vivo cell migration can substantially deviate from migration on 2D unconfined substrates (Davidson et al. 2014;Balzer et al. 2012).Indeed, during motion through tissues, ECM barriers, capillaries and lymph nodes, cells experience varying degrees of physical confinement and cell migration can thus be achieved with very different mechanisms (Lämmermann et al. 2008;Balzer et al. 2012;Bergert et al. 2015;Even-Ram and Yamada 2005).In particular, it has been observed that cell migration in 3D environments can occur even in the absence of focal adhesions, suggesting that additional mechanisms for adhesion and migration are possible (Lämmermann et al. 2008;Reversat et al. 2020).
Such adhesion-independent migration has been observed in 3D confined environments (Lämmermann et al. 2008;Friedl and Bröcker 2000;Friedl et al. 2001;Fraley et al. 2010;O'Neill et al. 2018), using different cell lines and technologies.For cells of different types (dendritic cells in Lämmermann et al. (2008); leukocytes in Reversat et al. (2020); breast carcinoma, pancreatic carcinoma, and human osteosarcoma cells in Balzer et al. (2012)), it has been observed that migration in 3D environments may not require myosin-mediated contraction and that inhibitors of integrins do not hamper migration through channels leading to cell confinement, although these treatments can hinder and even prevent motility in wider channels leading to unconfined migration.Considering leukocyte migration (Reversat et al. 2020), on the one hand, it was shown that leukocytes do not migrate when confined between two parallel flat plates in the absence of adhesion.On the other hand, leukocyte adhesion-free motion is possible when supporting pillars or microfabricated structured channels are placed between these plates, under the conditions that the pillar size and spacing-or the characteristic length of the sidewall structure-match roughly that of the cell length (Reversat et al. 2020).No cell migration is observed when the cell is confined between flat parallel plates in two directions, using the experimental set-up and the cell line reported in Reversat et al. (2020).
Even though the origin and transmission of propelling forces during focal adhesionfree migration are not fully understood, all these findings (Lämmermann et al. 2008;Balzer et al. 2012;Reversat et al. 2020) indicate that 2D is essentially integrindependent and that adhesion-free motility relies on a structured physical confinement, only achievable in a 3D setting, that can induce cytoskeletal alterations reducing the dependence of cell motion on the adhesion-contraction force coupling.In the absence of adhesions, non-specific transient interactions between transmembrane proteins and the substrate could generate friction that converts protrusive actin cortex flow into cell movement (Bergert et al. 2015;Hawkins et al. 2011;Lämmermann et al. 2008).It has also been observed that confined migration depends largely on microtubule (MT) dynamics and might persist even when F-actin is disrupted (Balzer et al. 2012;Stroka et al. 2014;Li and Sun 2018).
Although increasing levels of confinement can trigger transitions from integrinbased towards adhesion-independent migration modes in many cell types (Friedl et al. 2001), in the absence of matrix protease production, a too strong confinement either decreases or even prevents migration, due to cell stiffness (Wolf et al. 2003(Wolf et al. , 2013)).In particular, while the cytoplasm is very flexible and the cytoskeleton can actively remodel to undergo large deformations and penetrate small openings, the cell nucleus is normally 2-10 times stiffer than the surrounding cytoplasm and, with a typical diameter of 3-10 µm, occupies a large fraction of the cellular volume and is usually larger than many of the pores encountered in the extracellular environment (Davidson et al. 2014;Wolf et al. 2013;Cao et al. 2016).Thus, the nucleus should undergo substantial deformations when the cell moves through 3D constrictions, and it may constitute a rate-limiting factor during non-proteolytic migration of cells (Davidson et al. 2014;Wolf et al. 2013).
To understand the bio-physical and mechanical factors involved in the process of cell migration, many mathematical models have been proposed in the past decades (Jilkine and Edelstein-Keshet 2011;Holmes and Edelstein-Keshet 2012;Dreher et al. 2014;Danuser et al. 2013).Specifically, there have been abundant works related to cell migration on 2D substrates, either modelling the membrane mechanics and its signalling activity (Elliott et al. 2012;Hecht et al. 2011) or describing in detail the cytosol dynamics (Shao et al. 2010;Recho et al. 2013Recho et al. , 2015;;Manhart et al. 2015).However, coupled models, including the cytosolic machinery and membrane dynamics, have received little attention, even though they are critical to understand cell migration (Dreher et al. 2014;Danuser et al. 2013;Giverso and Preziosi 2018;Moure and Gomez 2018).Furthermore, most of these models have focused on 2D adhesiondependent cell motility, in which cells extend a stationary lamellipodium at the leading edge.Even in models accounting for amoeboid motion and which have been extended to model 3D confined migration (Moure and Gomez 2016, 2017, 2018), the cell motility substantially relies on adhesion, on acto-myosin protrusion-contraction, and on cell capability to sense an external field through membrane receptors.On the contrary, adhesion-independent migration inside constrained 3D environments has received less attention and mathematical models have started to tackle this interesting mechanism only recently.In particular, a simplified two-dimensional model for focal adhesion-independent cell swimming, based on the flow-friction driven force transmission, has been proposed by Wu et al. (2018) and Stotsky and Othmer (2022), while in Kaoui et al. (2008) the motion of closed phospholipid membranes suspended in a nonlinear shear gradient of a plane Poiseuille flow was investigated numerically in two dimensions.A possible explanation of the chemical signalling activity regulating adhesion-independent migration has been advanced in Elliott et al. (2012), using a system of reaction-diffusion equations and assuming a Turing instability to model a polymerization pattern on the cell surface, which drives the formation of pseudopods.We remark that the models for cell motion that introduce a friction coefficient between the cytoskeleton flow and the substrate to represent adhesion (Tawada and Sekimoto 1991;Giverso and Preziosi 2018;Farutin et al. 2019;Chelly and Recho 2022;Loisy et al. 2019) could as well be used to describe non-specific sliding friction (Farutin et al. 2019).Indeed, even though conceived for modelling the specific integrin-based adhesion, they can be adapted to describe transient interactions between the membrane and the substrate or the surrounding fluid, by performing appropriate calibration of the friction term.Recently, a couple of mathematical models have also investigated the non trivial limit of a vanishing friction coefficient with respect to other internal dissipative processes (Chelly and Recho 2022;Loisy et al. 2019;Le Goff et al. 2020), demonstrating that motility can still occur and making the models non-specific to the adhesion properties of the cell with its environment.However, since most of these models are interested in determining the minimal ingredient for the onset of cell motion, they are solved in a 1D setting (Farutin et al. 2019;Loisy et al. 2019;Le Goff et al. 2020).Furthermore, in all these cases (Kaoui et al. 2008;Wu et al. 2018;Stotsky and Othmer 2022;Elliott et al. 2012;Moure andGomez 2017, 2018;Chelly and Recho 2022;Loisy et al. 2019;Le Goff et al. 2020), the presence of the nucleus as a limiting factor for cell migration and the effect of confinement are not taken into account.
The influence of nuclear deformations on the whole process of cell migration was included in some recent works.In Moure and Gomez (2020) the role of the cell nucleus was studied using a computational model of a fish keratocyte, but the model is specifically conceived for 2D cell migration and thus it cannot be used for 3D confined migration.On the other hand, Cao et al. (2016) develop a chemo-mechanical model to study the nuclear strains and shapes, its plastic deformation and the threshold for the rupture of its envelope during migration through confined interstitial spaces.In Lee et al. (2017), a 2D model for cell migration through a dense network of host cells was proposed to reproduce glioma cell invasion.The moving cell is represented by two elastic closed curves, an inner curve corresponding to the nucleus of the cell and an outer curve corresponding to the cell basal membrane, whereas non-moving cells are represented by a single elastic curve.In Chen et al. (2018), the deformations of the cell and nucleus during invasion through a dense microenvironment were simulated incorporating stochastic processes and uncertainties in the input variables were evaluated using Monte Carlo uncertainty quantification simulations.These models Lee et al. (2017), Chen et al. (2018) and Cao et al. (2016) were able to reproduce correctly the hourglass cell and nucleus deformation observed in biological experiments, by relying on an external chemical factor.
In this work, we build on top of Jankowiak et al. (2020) to develop a simplified framework to study whether adhesion-free migration could be driven by simple mechanical features.We also test the influence of cell nucleus mechanical properties in the determination of the physical limit of cell migration.Thus, we propose a model of force generation during adhesion-independent cell migration in confined environments, taking into account the flow-friction driven force transmission, the cell membrane polymerization and the nuclear deformations.The cell is modelled by two membranes, an outer one representing the cell membrane and an inner one representing the nucleus.The two membranes are connected by microtubules, responsible for the nucleus location inside the cell.The renewal of the actin network underneath the cell membrane is modelled by the evolution of the mass distribution along the membrane, with a source (resp.sink) term at the front (resp.back), while also taking into account the conservation of the centre of mass.The model is used to simulate cell motion inside channels with structured walls with wavelengths ranging in the order of magnitude of the cell and nucleus diameters.We present the mathematical model in Sect. 2 and the numerical scheme we used to solve the systems of equations in Sect.3. Finally, in Sect.4, we present and discuss the numerical results.

The Mathematical Model
In this section, we present the continuum model ingredients for adhesion-free cell migration in domains containing rigid obstacles with a given geometry.Motivated by the experimental setup of Reversat et al. (2020), where the cell is confined between flat top and bottom surfaces and structured side walls, we choose a two-dimensional model, representing the projection of the three-dimensional set-up along the vertical directions, i.e. the one perpendicular to the flat walls.We consider the cell composed by two main compartments, the cytoplasm and the nucleus, both surrounded by membranes.The nuclear and cellular membranes can be represented as closed curves of R 2 .We identify the cell cortex with the cell membrane, and describe their ensemble as a single curve.This means that we do not model detachment and reattachment events between the cortex and the membrane, which are known to occur in certain cells.The cell cortex is schematically represented by a lipid bilayer and a complex underlying network of actin filaments.It is assumed to be elastic and subjected to a pressure differential force acting in its outward normal direction.The renewal of the actin network composing the cellular cortex is a key ingredient during many cellular behaviours and is here modelled by deposit and removal of material along the cortex.In rough biological terms, this corresponds to an imbalance between polymerization and depolymerization of some parts of the actin filaments and gives the cell a preferential direction of movement.The actin is then transported in the cell to a new location where it polymerizes again (actin treadmilling).This transport mechanism needs to be taken into account to ensure conservation of momentum in the absence of external forces.In our model, this is done by including an additional reaction force, as detailed in Sect.2.2.All these mechanical contributions on the cortex balance with frictional effects from the surrounding fluid (both external and internal to the cell), which are not explicitly described.
Finally, the nuclear membrane can be thought of as a double phospholipid bilayer with an associated mesh of intermediate filaments forming the nuclear lamina which stabilises the nucleus and provides some resistance to bending and tension.The nuclear envelope is eventually subjected to the differential pressure between the cytosol and the interior of the nucleus.Also in this case, we represent the ensemble of the nuclear membrane and lamina with a single curve.We use either the term nuclear envelope or nuclear membrane to refer to the whole structure of membranes and lamins surrounding the nucleus, as well as either the term cell membrane or cell cortex denote the outer membrane with associated actin cortex.

Description of the Model and Notation
Since we consider 2D-projections along the directions perpendicular to the flat plates of the 3D set-up used in Reversat et al. (2020), we set ourselves in R 2 and consider the cell cortex and the nuclear membrane represented by the time-dependent closed curves (t) and n (t) respectively, where T > 0 is some fixed maximal time and M > 0 is the fixed total amount, or mass, of actin along the membrane.The space variable s (resp.σ ) belongs to T M = R/MZ, which can be thought of as the interval [0, M] where 0 and M are identified, so that, for any fixed t ∈ [0, T ], the continuity of X (t, •) and Y (t, •) is enough to enforce the closedness of the curves (t) and n (t).Note that X and Y are not assumed to be arc length parameterizations.The time derivatives are denoted by ∂ t X and ∂ t Y whereas the space derivatives (along the curves) are ∂ s X and ∂ σ Y , respectively.In some sense, the variable s counts the mass of actin along the curve and therefore tracks Lagrangian particles.More precisely, for any non empty interval [s 1 , s 2 ] with 0 ≤ s 1 < s 2 ≤ M, the amount of actin on the corresponding piece of cortex, i.e. between X (t, s 1 ) and X (t, s 2 ), is |s 2 − s 1 |.Therefore, X encodes not only geometric information (the shape of the curve), but also information about the distribution of actin inside the cortex.By writing the problem in terms of s (as opposed to arc length), we can describe the time evolution of both these quantities with a single equation for X .Writing the equation in terms of the arc length ∈ [0, L] makes it harder to deal with a time-dependent length L, as is the case here, especially for the numerics.Therefore, we work with the time independent T M s.For the sake of completeness, let us highlight the link between the variable s, the arc length and the actin density on the cortex, which we will denote by ρ.As is standard, the length between two points X (t, s 1 ) and X (t, s 2 ) is given by As already mentioned, the mass of actin between X (t, s 1 ) and If we now consider s to be a function of the arc length , we have, with a slight abuse of notation, that . By a change of variables, this is compatible with In the following, we take M = 1, implying that the actin mass is normalised with respect to a reference total mass, so that s ∈ T 1 .
Assuming that the curves are smooth enough, we denote by τ (t, s) and n(t, s) the unit tangent and unit outward normal vectors to the cell membrane curve at X (t, s) and with T (t, σ ) and N (t, σ ) the unit tangent and unit outward normal vectors to the nuclear membrane curve at Y (t, σ ).Assuming positive orientation of the parametrization, we have with the convention (a, b) ⊥ = (−b, a).A sketch of the notations employed in the paper is illustrated in Fig. 1.Finally, let (t) ⊂ R 2 be the set bounded by (t) and n (t) ⊂ R 2 the set bounded by n (t), i.e. (t) = ∂ (t) and n (t) = ∂ n (t): for biological consistency we have to guarantee that the nucleus is located inside the cell, i.e. n ⊂ .During cell motion and other biological processes (e.g.development, mitosis, fertilisation, …), Fig. 1 The parameterization and associated vector quantities, along with the representation of the microtubule network geometry, with the centrosome X c linked to the centroid of the nucleus Ȳ .The pink region shows where microtubules are present, the corresponding anchoring points on the membrane are represented in bold stroke.The bottom arrows indicate the orientation of the curves.The elastic link between the centrosome and the nucleus (coloured in cyan) is also represented, as well as the microtubular force F MT , both of which will be detailed later on (Color figure online) the positioning of the nucleus is of paramount importance in the establishment of cellular architecture (Tran et al. 2001).Nuclear positioning is generally dependent on some cytoskeleton constituents, mainly microtubules (MTs), which are dynamic polymers of tubulin, and intermediate filaments, composed of a family of related proteins having common sequence and structural features.Microtubules originate from the MT-organising centre (MTOC), with functions that include microtubule nucleation, stabilisation, and/or anchoring.The best-studied MTOC is the centrosome, which is found in many animal cells.The centrosome is connected to the cell nucleus through MTs, and associated intermediate filaments forming a ring network around the nuclear envelope, as well as to the cell actin cortex, via MTs.Thus, the MT structure and its related centrosome couple the nucleus to the cellular envelope, and play a fundamental role in providing structure and shape to cells, in determining cell migration direction and persistence, and in locating the cell nucleus (Fruleux and Hawkins 2016;Gundersen and Worman 2013;Beadle et al. 2008;Friedl et al. 2011).Microtubules constantly switch between growing and shrinking states, through assembly and disassembly of tubulin monomers at their ends, in a process termed dynamic instability.They are able to generate pushing forces over the cell membrane during their assembly process, and pulling forces during their disassembly process (Laan et al. 2008).In this work, we disregarded the description of the dynamic growth/shrinkage of the MTs through the addition and detachment of monomers, but the MT structure is built at each time t, depending on the position of the centrosome and the cell cortex.In this scenario, without describing the incremental growth of each filament over time, the MT structure evolves and it is not fixed once and for all.By adopting a continuous setting, at each instant of time, the MTs are assumed to be homogeneously distributed around the centrosome, located in X c (t), to all points on the cortex that can be connected to the centrosome by a line segment, lying inside the cell.The region of the cell in which MTs can be defined is coloured in pink in Fig. 1.
With that in mind, it is possible to define, for every time t, the microtubule anchoring points on the membrane cortex as the first intersection of the half line starting at X c (t) with angle θ ∈ [0, 2π) with the cell cortex.Formally, we define the map where λ MT (t, θ) := min {λ ≥ 0 : X c (t) + λe θ ∈ } and e θ = cos(θ ) sin(θ ) .
This map MT is well defined if X c (t) ∈ (t) and is surjective if (t) is star-shaped with respect to X c (t).Then, the cortex anchoring point of the microtubule located at an angle θ is defined as In Fig. 1, the portion of the cortex on which MTs can anchor is highlighted in bold stroke.Note that, with a slight abuse of notation, MT can also be seen as a map onto Then, the segment [X c (t), X (t, s MT (θ ))] represents a microtubule at a given instant of time.The construction of some representative MTs is illustrated by the purple segments in Fig. 1.We remark that MTs can be drawn in the whole region of the cell coloured in pink in Fig. 1 defining a continuous MT structure.Therefore the MTs structure homogeneously spans all the angles around the centrosome, and the distribution of the length of MTs can vary over time depending on the location of the centrosome and the cell membrane points.Since the 2D representation used in this paper stands for a projection of the 3D cell along the vertical direction, we assume that the MTs and the centrosome can be built also in the nucleus region (cyan area in Fig. 1), representing filaments extending either underneath or above the nucleus.Finally, we observe that in this description we disregarded the deformation of MTs and we cannot capture the MTs' bending at the cell membrane (Geisterfer et al. 2020).However, the description of such phenomena will require a deeper mechanical characterisation of the MTs and will call for a 3D model in order to correctly represent the geometry of the MT filaments.
In the following, under the setting depicted above, we will derive the equations describing the evolution of the cell cortex (Sect.2.2), the MTs structure (Sect.2.3), and the nucleus membrane (Sect.2.4).The dependence of the different dependent variable functions on the independent variables in their arguments will be omitted whenever possible.

Evolution of the Cell Membrane
Concerning the evolution of the cell membrane, we refer to the model proposed in Jankowiak et al. (2020), properly modified in order to take into account the presence of the nucleus and the MT structure.The evolution of the actin density, describing the active component of the model due to the heterogeneity of the polymerization rate across the cortex, is given by Jankowiak et al. (2020): where f (t, ) is the rate of actin density increase ( f (t, ) > 0) or decrease ( f (t, ) < 0).As done in Jankowiak et al. (2020), we assume that the total amount of actin in the cortex is kept constant and that the cell polarization and subsequent actin polymerization manifests itself by a local imbalance producing a net increase of actin density close to the cell front, balanced by a decrease close to the rear of the cell.Since it is assumed that the total amount of actin in the cortex does not change in time, we require The mass transfer rate, or polymerization rate, r pol ≥ 0 is then defined as where ( • ) ± denotes the positive and negative parts, respectively.
In practice, we assume that the cell is polarized in a given direction e p ∈ S 1 , and that for each time t, there are unique s back (t), s front (t) ∈ T 1 (see the left of Fig. 2) so that A reasonable choice for f is then a function with its (non negative) maximum at X (t, s front (t)) and (non positive) minimum at X (t, s back (t)).
For what concerns the MT endpoints density on the cell cortex, we can define To derive Eq. ( 7), we exploit the fact that Then, the evolution of the cell cortex is determined by the Newton's Second Law, written in an overdamped regime, i.e. the friction of the surrounding fluid balances with all other contributions, leading to the following balance law for X : where the friction coefficient in front of the time derivative is set equal to unity, without loss of generality, by an appropriate choice of the time scale, as done in Jankowiak et al. (2020).We remark that the description of the motion of the fluid both inside and outside the cell is here neglected.This assumption seems reasonable, at least referring to the biological experiments performed by Reversat et al. (2020), since we study the motion of a cell confined inside a channel where the hydrodynamic interactions are not the leading cause of cell motion.In Eq. ( 8), the velocity of the cell cortex relative to the laboratory coordinates is given by the material derivative D X(t, s)/Dt and not by the partial derivative ∂ X (t, s)/∂t.This is because of the implicit dependence of s on t which is introduced by the time evolution of the density ρ in (5).In Fig. 2, we sketch why, for a nonzero polymerization rate r pol and for fixed s, X (t, s) does not track a material point.The material derivative is then defined as where we used (1) and ( 5).
The second term on the l.h.s. and the first term in the r.h.s. of eq. ( 8) represent, respectively, the friction force generated by microtubule-binding complexes sliding on the cortex and the in-line force due to MTs elongation (as they will be detailed in Sect.2.3).Both terms act on the portion of the cell membrane where MT anchor points can be defined and therefore they are weighted by ρ MT , defined by Eq. ( 7).The second and third terms on the r.h.s. of Eq. ( 8) take into account, respectively, the contact between the cell cortex and the nucleus (see Sect. 2.4 for a detailed description of F cont ) and the one between the cell and the channel wall.In the following, we assume that F wall derives from a potential, a possible choice of which is proposed in Sect. 4. The compensating force F comp in the fourth term on the r.h.s. of Eq. ( 8) must be chosen (for each time t) so that the centre of mass is fixed if the effects of the confinement, nucleus and microtubule structure are removed and if the motion of actin along the cortex is perfectly balanced by the treadmilling, without any internal dissipation (see Jankowiak et al. (2020) for more details on the derivation).In this show how the material points at t and s = 0 and s = s 1 are mapped at t + ε.Because of the polymerization, ∂s/∂t > 0 in the upper part of the curve between the two yellow regions, so that, at t + ε, the material point which was located at case, we obtain after integration where F T (s) is defined through Eq. ( 9) and it is related to the actin transport due to polymerization and depolymerization.Finally, the last term on the r.h.s. of Eq. ( 8) takes into account the forces acting on the cell due to the cortex mechanical behaviour and the pressure difference in and out of the cell and it can be obtained from the cell membrane energy The membrane energy E c comprises an elastic term representing the cell and cell cortex elasticity E el and a term related to the existence of a differential pressure E p .More precisely, the associated elastic energy E el is composed of two parts, the first one E (1)  el is related to the response of the cell membrane to tension and the second one E (2)   el is related to the response of the cell to deviations from its target area.The introduction of an elastic constraint on the cell membrane length and cell area is in agreement with prototypical models of membrane and cortex elasticity using elastic springs linking different parts of the cell (Barnhart et al. 2010;Du et al. 2012;Recho and Truskinovsky 2013;Kuchnir Fygenson et al. 1997).We observe that the elastic response of a cell may be associated with both the actin cortex and the phospholipid cell membrane.Since the actin cortex undergoes a constant renewal over a timescale of 30-100 s (Rubinstein et al. 2009) and bulk elastic stresses inside the cortex are relaxed over a time-scale of 1-10 s (Rubinstein et al. 2009;Mofrad 2009;Recho et al. 2015), which are much shorter than the characteristic timescale of motility experiments, the elastic behaviour of the outer curve is mainly associated with the mechanical response of the cell membrane itself and not of the actin cortex.
Formally, E (1) el reads where k c is the mechanical parameter representing the cell membrane stretchability.This choice is motivated by some reasoning at the discrete scale.Indeed, we consider the following very simple model of the cortex: we think of the cortex as a (closed) chain of n individual point masses (X i ) 1≤i≤n , each of mass M/n.They are linked by Hookean springs of stiffness kc and of equilibrium length 0 , so that the potential energy of the spring between X i and X i+1 has the expression The total potential energy is obtained by summation over i: By considering the scaling kc = k c −2 0 with 0 ∝ n −1 , we can take the (formal) limit as n → ∞ to recover E (1) el .Since the index i counts the "mass", so does the continuous variable s, as explained at the beginning of Sect.2.1, it is possible to obtain Eq. ( 11).One can also write E (1) el in terms of the arc length: so that the local energy density is convex, with its minimum in ρ ≡ 1, and becomes large for small or large values of ρ from the reference density, i.e. either ρ 1 or ρ 1.We remark that, even though displacement does not appear explicitly, the contribution of E (1)  el tends to keep |∂ s X | (or equivalently, ρ), and thus the length L, close to one, thanks to Eq. ( 2).If one assumes that the actin cortex is continuously laid out (i.e.polymerized) with a density ρ = 1, ρ −1 can then be interpreted as some form of displacement, which explains the form of the right-hand side in the equation above.The factor ρ can be understood as follows: for a given displacement, the local energy density grows linearly with the density of actin filaments.The above expression is similar to that of the standard potential energy for an elastic ring, except for the factor |∂ s X | −1 = ρ.In other words, the local energy density is also proportional to the actin density ρ, which is consistent with the discrete model: for fixed n, as the point masses X i get far from one another, the number of springs in a given length interval decreases, and so does the force they exert.Furthermore, we include an elastic constraint on the cell area | | (which would correspond to the volume in three-dimensions).In particular, according to previous works (e.g.Kuchnir Fygenson et al. 1997), we assume that the elastic energy of the cell is minimised if the cell area is equal to a given target area A * c , i.e.
where μ c is the mechanical parameter representing the elastic resistance of the cell bulk to variations of its area, given by the measure of the domain , | |.Physically, this models the resistance to compression of the organelles in the cytoplasm.
Concerning the term E p , because of osmotic effects, we suppose that the cell is subject to an internal cytoplasmic pressure, which results in a force in the direction of the normal to the curve.The force intensity per unit length is assumed to be uniform in space and constant in time, so that the associated energy E p is where p > 0 is the constant excess of pressure inside the cell, with respect to the extracellular pressure.

Evolution of the Microtubules' Structure and the Centrosome
As depicted in Sect.2.1, microtubules (MTs) and the associated centrosome are the principal coupling mechanisms between the cortex and the nucleus.Microtubules are known to generate forces to position and shape the cellular organelles, and in particular the nucleus (Mofrad 2009).A number of observations (Soheilypour et al. 2015;Mofrad 2009;Stamenović et al. 2002) suggest that among all cytoskeletal components, MTs play a critical role in carrying compressive loads, behaving as passive compressionsupporting elements that maintain cell shape.However, when MTs are anchored to the cell cortex, through dynein motor proteins, they are able to generate pushing forces over the cell membrane during their assembly process, and pulling forces during their disassembly process (Laan et al. 2008).Since the positioning of the nucleus and the centrosome has been found to principally depend on MTs pushing/pulling forces (Laan et al. 2008), we are here interested in giving a simplified mathematical description of such a force, directed along the microtubule axis.
In the mathematical setting put forward in this work, at every instant of time, a continuous structure of MTs, homogeneously distributed around the centrosome, located in X c (t), connects the centrosome to the points on the cortex s MT (t, θ) = MT (t, θ), ∀θ ∈ [0, 2π) where MTs can anchor (see the bold line in Fig. 1).
The microtubule push/pulling force F MT (t, θ) is, thus, directed along the line segment [X c (t), X (t, s MT (θ ))], that represents the microtubule (see the reference purple segments in Fig. 1).Many works in the literature focus on MT mechanical response to bending (Mizushima-Sugano et al. 1983) and radial indentation (Schaap et al. 2006), but less is known on their response to elongation.However, some works highlight interesting behaviour of the MTs with respect to forces directed along the MTs major axis (see Laan et al. 2008;Soheilypour et al. 2015), which is responsible for the positioning of the nucleus and the centrosome.In particular, in Laan et al. (2008) it was assumed that the in-line force exerted by a MT is related to its polymerizing activity, whereas in Soheilypour et al. (2015) the buckling of long MT bundles was observed, which can strongly reduce the exerted tension.In this work, we consider that the modulus of the MT force, directed along the MT major axis, is a function of the microtubule length |X c (t) − X (t, s MT (θ ))|, to be specified, i.e.
were e θ is the outward unit vector representing the direction of a MT.In the following, we will consider , so that smaller MTs, that can polymerize the most, exert the higher forces, whereas the force decreases with the MT length, due to buckling instability.We observe that, since f * MT > 0, the force acting on the MT is always directed along −e θ , whereas the force exerted by each MT on the cell cortex is pushing against it, being directed along e θ .We remark that the microtubule force acts on the portion of the cell cortex where ρ MT = 0 [see Eq. ( 8)] and on the centrosome, where MTs originate.The total resultant of all microtubule forces acting on the centrosome is then, In addition to the in-line forces due to elongation, the cortex endpoints of the microtubules are also subject to a friction force, caused by the sliding of the binding complexes on the cortex.This force is directed opposite to the velocity of the centrosome relative to the cortex.The friction also generates a torque, with angular velocity ω.The motion of the MT structure with respect to the angular velocity can be regarded as rigid in the sense that all MTs rotate with the same angular velocity.However, since the MT structure does not simply rotate but also resizes according to the cell shape and the microtubule organising centre location, it is not a standard rigid body rotation.As already pointed out, this is of course an approximation of the far more complex biological reality.
Defining v(∂ t X c , ω, θ) the speed of the centrosome relative to the cortex endpoint of the microtubule of angle θ , we have The first term on the r.h.s of Eq. ( 15) represents the velocity of the centrosome, whereas the second term, between square brackets, stands for the velocity of the cortex point where the microtubule is attached.Since the MT structure is redrawn in the cell at each instant of time, without explicitly tracking the motion, elongation and/or breakage of each single MT, in order to properly define the velocity of the MT on the cortex, we have to consider the anchor point of the MT onto the cell cortex whose position is well defined and can be differentiated in time, instead of the end of the MT, which is not explicitly tracked.We remark that we consider ∂ t X (t, s MT (t, θ)) as opposed to the total derivative D Dt X (t, s MT (t, θ)), since we want to describe the relative velocity of the centrosome with respect to the cell membrane and not with respect to the laboratory frame.Furthermore, the third term on the r.h.s. of Eq. ( 15) corresponds to the sliding due to the rotation of the MT structure at each instant of time, which is obtained considering the relation between the linear and the angular velocity.
For the benefit of subsequent computations, the relative velocity v can be also decomposed along the directions τ and e θ : where τ and e θ are the corresponding oblique projections: and Locally on the cortex, the friction of the microtubule against the cortex generates a force where k τ and k e θ are the friction coefficients in the corresponding directions.In the following, for the sake of simplicity, we will take k e θ = k τ , but the model can be easily generalised to the case k e θ = k τ .By writing F int := FX c + FMT the sum of forces acting on the microtubule structure, which include FMT (the forces directed alongs the MTs), given by Eq. ( 14), and FX c , the forces acting on the centrosome itself (that will be detailed in Sect.2.4, see Eq. ( 27)), we have the following force balance in R 2 : where we used the shorthand X θ = X (s MT (θ )).
Since the microtubule structure has zero moment of inertia w.r.t.X c , we also have the following balance of torques which can be rewritten by decomposing v on e θ and e ⊥ θ : Gathering ( 16) and ( 17), we get the system where where I 2 is the identity matrix in 2D.We remark that, by Jensen's inequality 1 , we have The equality occurs only if X θ reduces to a single point, so we can assume that the determinant is always positive.Then, system (18) can be solved to get the angular velocity ω of the MT structure and the velocity of the centrosome: 1 Jensen's inequality states that φ

Evolution of the Nuclear Membrane
We consider the nucleus material to be harder to deform than the rest of the cell and with a negligible relaxation time.The nuclear membrane has a certain mechanical behaviour and, as previously stated, it is connected to the centrosome, which contributes to the positioning of the nucleus inside the cell.Furthermore, the nucleus cannot cross the cellular membrane and therefore interacts with it through a contact force F cont,n .Then, calling F n the forces internal to the nuclear membrane, F X c N the force acting on each point of the nuclear membrane due the interaction between the nucleus and the centrosome, F cont,n (σ ) the contact force with the cortex, the Newton's Second Law for the nuclear membrane, neglecting inertial terms, reads The meaning and derivation of the different terms on the l.h.s of Eq. ( 21) is depicted below, whereas the term on the r.h.s.describe the friction force acting on the nucleus.
In the following, we will set the coefficient h v equal to unity, assuming that the viscous coefficient representing the interactions between the nuclear membrane and the cytosol is the same as the one representing the interaction between the cell cortex and the intra-and extra-cellular fluid.This assumption seems reasonable when the cell migrate in an adhesion-independent manner, whilst in adhesion-mediated migration the friction term in the cortex equation should also take into account the formation and breakages of adhesion points with the surrounding environment and it is certainly different from the viscous coefficient representing the interactions between a cellular/nuclear membrane and a fluid.For what concerns the l.h.s of Eq. ( 21), the force F n related to the mechanical behaviour of the membrane can be derived from the nuclear membrane energy, E n , through the relation where the energy E n is the sum of all energies related to the constitutive mechanical behaviour of the nuclear envelope, i.e.
where d = |∂ σ Y |dσ is the arc length element and d A the area element.The first term in Eq. ( 22) represents the energy associated to the membrane bending, k b is the bending modulus of the nuclear membrane, K the local curvature and K 0 the characteristic (or spontaneous) curvature, which is assumed to be zero in what follows, according to Kaoui et al. (2008).The second term represents the tensile stress acting on the membrane and λ can be thought as the nucleus surface tension.p n is the difference between the pressure in the cytosol and the pressure inside the nucleus.Finally, the last term represents the volumetric elastic constraint associated to changes in nucleus area A n relative to a defined target area A * n and it has the same form as the constraint introduced for the whole cell in Eq. ( 12).
Taking all these contributions into account, the computation of the variation δ E n of the nuclear membrane energy function (see Section A) leads to the following nuclear membrane mechanical force where N is the outward normal to the nucleus (see Fig. 1).Then, to derive a mathematical expression for the second term on the r.h.s of Eq. ( 21), we assume that the centrosome is linked to the nuclear membrane points thanks to cytoskeletal filaments, namely microtubules and intermediate filaments (Cooper 2000).Intermediate filaments form an intricate ring network that surrounds the nucleus of most cells and extends in the cytoplasm, where they associate with the other elements of the cytoskeleton, such as microtubules (Cooper 2000).Intermediate filaments attached to the nuclear envelope, along with microtubules, serve to position and anchor the nucleus within the cell, to redistribute forces acting on the nucleus and to provide a scaffold that integrates the components of the cytoskeleton (Cooper 2000).Therefore, assuming that each of these cytoskeletal coupled structures behaves as a sort of spring and that the network formed by intermediate filaments surrounding the nucleus leads to a uniform distribution, along the length of the nuclear membrane, of the forces coupling the centrosome to the nuclear envelope, we have where k e is the elasticity of each virtual link between the centrosome and the nuclear membrane point, and L n is the length of the nuclear membrane at time t, which is given by Since the 2D setting that we are considering is to be interpreted as the projection of the 3D geometry, we assume that the rest length for the connection between the nucleus centroid and the centrosome is equal to zero.This condition represents the 3D situation in which the centrosome is on the top or bottom of the nucleus.It is placed at the centre of the 2D projection of the nucleus.If one assumes that k e is constant, it is possible to rewrite (23) as where Ȳ is the centroid of the nucleus (see Fig. 1) defined as Under the same assumption, it is also possible to define the total force FX c acting on the centrosome due to all centrosome-nucleus interactions, which is thus part of F int in Eq. ( 16): We observe that the presence of an elastic coupling between the nucleus centroid and the centrosome, specified in Eqs. ( 25) and ( 27), is in agreement with the datadriven theoretical approach proposed in Brückner et al. (2022).Specifically, Brückner et al. ( 2022) develop a model for protrusion and polarity dynamics in confined cell migration, combining experimental data inference with a mechanistic approach, based on an elastic coupling between the cell nucleus centroid and the center of the cytoplasm protrusion, combined with a non-specific friction acting both on the cell and the nucleus and a proper description of cell polarity.Finally, as previously mentioned, the nucleus is constrained to remain within the cell, since it cannot cross the cortex.This is guaranteed by including a penalization, represented by the last term on the l.h.s of Eq. ( 21).The contact force is given by where W cont is a decreasing function with compact support.By symmetry, we specify the corresponding force acting on the cortex, F cont , appearing in (8).
We remark that the inclusion of the contact constraint is essential in order to prevent the penetration of the nucleus and the cell cortex, for any range of the model parameters and channel size.Indeed, the inclusion of the MTs force given by Eq. ( 13) is not enough, since it does not directly act on the nuclear membrane points, but its resultant (see Eq. ( 14)) is applied to the centrosome.The location of the centrosome, in turns, affects the positioning of the center of mass of the nucleus, but without any contact force between the nucleus and the cell cortex, the penetration between some regions of the nucleus and the cortex could take place.The specific functional form used for W cont will be discussed in Sect. 4.

Numerical Discretisation
The whole system composed by Eqs. ( 8)-( 18)-( 21) is discretised in space using finite differences, the resulting system of ODEs is then solved using split step time stepping scheme of order one.At each time step, the new position of the nuclear membrane is computed using the explicit scheme described below.Then, the position of the cortex and centrosome are updated using a semi-implicit scheme.

Cortex, Centrosome and Microtubules Structure
We assume the total mass of actin on the cortex to be normalised to 1 and we introduce N 1 ∈ N grid-points for the discretisation of s (corresponding to the cortex) such that s i = i s, s = 1/N 1 for i ∈ {0, N 1 − 1}.Given a time step t > 0, we denote t j = j t.In what follows, subscripts (resp.subscripts) correspond to space (resp.time), and we define X j i :=X (t j , s i ).For the sake of legibility, indices or exponents are omitted if possible.For integrals, we use the midpoint rule, meaning that and The cortex elastic force, the bulk elasticity and pressure forces are discretised as where the (polygonal) area A j is computed as 123 The interaction force with the wall simply becomes whereas the transport term in (9) reads Following (10), the compensating force F comp must satisfy In the present work, we are mostly interested in the impact of the nucleus on the dynamics, so we can choose F comp,i independent of i for simplicity.The computation of MT (or equivalently s MT (θ )) and ρ j MT ,i requires the construction of the visibility polygon of the cortex from the centrosome.We do not discuss the construction here and refer to Lee (1983) and Joe and Simpson (1987) for details.The quadrature formulae and corresponding expression for MT are detailed in "Appendix B".
It remains to deal with the contact force between the cell cortex and the nuclear membrane, which we simply take as the discretisation of (29): The Eq. ( 18) for the angular velocity ω of the MT structure and the velocity of the centrosome is discretised similarly in time and space, using the quadrature formulae reported in the Appendix (see section B) to compute the integral of the different quantities.
For the time iteration, we use an implicit Euler scheme, so that for both the cortex and the centrosome we obtain a system of the following form at each time step: , where the r i , r c , r ω on the r.h.s.correspond to the discretized terms detailed in this section, and the matrices M are the corresponding jacobian matrices.The left-most matrix above corresponds to the discretization of the time derivatives.Since the problem is formulated in term of the angular velocity ω, whose governing equation do not involve derivatives in time, the last element of the diagonal is zero.

Nuclear Membrane
The evolution of the nuclear membrane is derived rewriting the proposed equation ( 21) with the formulation and the efficient discretisation scheme proposed in Mikula and Ševčovič (2004) and Beneš et al. (2009) appropriately adapted to our setting.The method basically consists in splitting the velocity guiding the evolution of the nuclear membrane, ∂ t Y , into its normal component, β, and its tangential one, α, so that ∂ t Y = β N + αT .The normal velocity is the one determining the variations in the nuclear morphology, whereas the tangential component has no impact on the shape of the evolving curve, but it governs the distribution of the nodes along the nuclear envelope.The evolution of the nucleus is then given by the following system (see the Appendix, Sect.C.1 and C.2, for further details) where is the arc length, K , ν and η = log |∂ σ Y | are, respectively, the curvature, the tangential angle and the logarithm of the local length of the nucleus envelope n at a point Y ∈ n , whereas W (Y ( )) is the potential combining the contact interaction with the cell cortex, the elastic constraint with the centrosome and the nuclear membrane surface tension.In Eq. ( 31) we take where n n u( ) d denotes the average of u over the whole curve n , which makes (32) a non-local equation, whereas ζ > 0 is a given positive constant, included in order to avoid nodes to concentrates on points, which would lead to poor approximation and eventually inversion of ill-conditioned matrices.
To write the corresponding discretisation of (31), we uniformly discretised the fixed parameterization interval [0, 1] in N 2 subintervals, each of equal length h = 1/N 2 and indexed by i ∈ {0, N 2 − 1}.For time, we use the same discretisation introduced for the cell membrane, so that the point Y (ih, j t) is written Then, the system of equations ( 31)-( 32)-( 33) is solved for the discrete quantities α The algebraic system determining the evolution of the nuclear envelope is reported in the "Appendix C.2", along with some comments on the derivation of the discretised equations.Finally, the discretised system of equations has been solved implementing a numerical code with Julia 2 (Bezanson et al. 2017).

Setup
In this section, we present the numerical results obtained solving the discretised system of equations presented in Sect.3. In particular, we consider channels with structured side walls of the following form: where f width represents half of the mean width of the channel, f β < f width is the amplitude of the oscillation of the wall, and f ω 0 is the pulsation of the sinusoidal channel (see Fig. 3 for an illustration of these quantities).We remark that for f β = 0, one recovers the flat walls geometry.We then choose where g ξ (x) = −min(ξ x − 1, 0) 2 log(ξ x) is a smooth barrier function with lim x→0 + g ξ (x) = +∞ and g ξ (x) = 0 for x > ξ −1 .Analogously, for describing the contact force between the nucleus and the cell cortex in Eqs. ( 28)-( 29), we use where k cont , ξ cont > 0 are parameters.
Concerning the polymerization, we follow Sect.2.2 and choose f as a super-Gaussian: f = exp(−( x 2 2w ) P )/C(t), where C is a normalisation factor.We choose w = 0.5 and P = 3.
Finally, we will consider the case in which the cell is represented • only by the cell envelope, as done in Jankowiak et al. (2020), for comparison; • by the cell envelope and the cell nucleus, linked together by the microtubule structure, as explained in the Sect.2.
The initial condition for the cell cortex is chosen as the evenly spaced discretisation of a closed curve C 0 which matches the side walls of the channel-albeit with a smaller 2 julialang.org.width-in order to have an initial cell area equal to the target area A * c .More precisely, it is the union of the following 4 curves: )} where x max 0 is chosen so that the area enclosed by C 0 is A * c .The initial condition for the nucleus is the circle n,0 centred on (π/2 f ω 0 , 0) and such that n,0 + B ξ −1 (0) ⊂ wall , where B ξ −1 (0) is the open ball of radius ξ −1 and + denotes the Minkowski sum.This construction is illustrated in Fig. 3

(right).
The spatio-temporal evolution of the cell and nuclear envelope for some benchmark simulations are reported in Sect.4.2, whereas the effect of the different parameters of the model on cell ability to move and its velocity inside sinusoidal channels is investigated in Sect.4.3.

Results: To Move or Not to Move
In this section, we investigate the ability of the model to reproduce cell migration inside a channel and the spatio-temporal evolution of the cell and the nucleus shapes.
We first consider the case of a cell positioned inside a channel with flat walls ( f β = 0): in this case, independently on the widths of the channel, no migration is observed in the simulations, in spite of polarization and corresponding cortex flow (see Fig. 4).This is due to the total balance in the transport of actin in the cell.If one only considers polymerization and the retrograde flow of actin, one expects forward movement.However, actin which is depolymerized at the back has to be transported towards the front.This mechanism, which is not modelled explicitly but taken into This result confirms that adhesion-free motility relies on structured confinement and it is in agreement with the experiments performed on leukocytes where retrograde cortical flow can be observed without motion of the cell (Reversat et al. 2020) and with the model described in Jankowiak et al. (2020), where the nucleus is not considered.
We then describe the motion of the cell inside a channel with sinusoidal walls.In this configuration, the cell takes an hourglass configuration with nucleus deformation which is often observed in vivo, so that channels of this type can be argued to mimic the inter-and extra-cellular space in which cells naturally migrate.
In this case, as it is known from biology, the capability of the cell to migrate inside the channel and its speed of migration are highly influenced by the presence of the nucleus, which is the stiffest part in the cell and can therefore remain stuck in the rear of the cell, preventing cell motion.In the present model, the resistance of the nucleus to deformations is highly controlled by the bending modulus k b and by the elastic areachange constraint μ n .Therefore, in Fig. 5 (and in the videos in the Supplementary materials) we report the evolution of the cell and the nucleus shapes at the same instant of time, for some typical simulations, obtained varying the parameters k b and μ n .Namely, we consider (a) the migrating cell without the cell nucleus, as done in Jankowiak et al. (2020) (Fig. 5a); (b) the migrating cell with a low bending modulus of the nucleus (Fig. 5b, with k  The other dimensionless parameters chosen in the simulation are summarised in Table 1.We observe that, for the particular choice of parameters set in the simulations  1 reported in Fig. 5a-c', the cell is able to migrate inside the channel even when the nucleus is explicitly modelled.The cell and the nucleus size and shape change during the motion as a balance of the cell area and nucleus area penalisations, their mechanical properties, microtubules and polymerizing forces, and the contact with the channel wall.In particular, it is possible to see that the cell first protrudes the cytoplasm to fill the maximum of sinusoidal spaces ahead the cell nucleus.The maximum cytoplasm extension is controlled by the cell target area, the cell membrane target surface, and by the mechanical parameters of the cell membrane.Since these parameters are kept fixed in the simulations in Fig. 5, the cytoplasm extension inside the channel is comparable for the different cases.When the cytoplasm completely fills the new space, the nucleus is pulled by the microtutubule structure through the constriction in the channel and it forms a bleb in the front until the nucleus is pushed inside the second constriction in the sinusoidal channel.When the nucleus has a low bending stiffness k b (see Fig. 5b), it acquires an hourglass shape by passing through constrictions in the channel and the formation of nuclear protrusions is evident both at the cell front and at the cell rear.However, for higher values of bending stiffness (see Fig. 5c) the nucleus shrinks to pass through constrictions and the formation of blebs and the development of an hourglass shape are less pronounced.Moreover, if we increase the nucleus area-change constraint, μ n , by keeping the value of k b fixed, the nucleus cannot shrink to pass through the constrictions and intense nuclear deformations, with blebs both in the front and in the rear of the nucleus, are observed (see Fig. 5c').In all these cases, when the nucleus fills the new sinusoidal space ahead the cell, the cytoplasm protrudes in the next available sinusoidal space and the process is repeated cyclically, allowing the cell to move forward inside the channel.
We remark that, although the cell motion occurs inside a simplified geometry, the nucleus hourglass deformation and the formation of blebs in correspondence of small openings in the extracellular space is a characteristic observed during cell motion inside intricate ECM (Wolf et al. 2007(Wolf et al. , 2013;;Beadle et al. 2008).
Furthermore, the importance of including the description of the nucleus becomes clear from simulations reported in Fig. 5d, in order to account for those situations in which the cell can remain trapped inside the channel, because the nucleus cannot deform and squeeze through the small opening of the channel, thus preventing the cell motion.Indeed, in this last case, even though the cytoplasm protrudes inside the channel and the nucleus is pulled by microtubules, the energy required to deform and bend the nucleus is too high (due to the high value of k b ), so that the nucleus gets stuck in the rear of the cell (see Fig. 5d).This finding is in qualitative agreement with a number of experimental works, such as Wolf et al. (2007Wolf et al. ( , 2013)), Rolli et al. (2010) and Beadle et al. (2008), where the cell migratory capability is associated with nuclear deformations, and the existence of a critical ECM gap size below which cell migration is entirely hampered in the absence during non-proteolytic has been observed.Such a critical size was termed "the physical limit of cell migration" (Wolf et al. 2013).In particular is has been observed that, when passing through constrictions, the nucleus shape can strongly deviate from the spherical one.Specifically, the shape of the deformed nucleus inside regular cylindrical channels can be approximated either by a prolate ellipsoid (Versaevel et al. 2012;Friedl et al. 2011) or by a cigar-like shape (Friedl et al. 2011), whereas inside structured channels or complex 3D extracellular matrix, the shape of the nucleus can highly vary, from a regular hourglass shape to more irregular nuclear conformations (see the experimental pictures reported in Friedl et al. 2011;Davidson et al. 2020Davidson et al. , 2014)).
Moreover, the numerical results in our work are in agreement with other mathematical models, dealing with the influence of the mechanical properties of the nucleus on the cell's ability to migrate in channels composed of extracellular matrix (Giverso et al. 2014(Giverso et al. , 2018;;Scianna et al. 2013;Scianna andPreziosi 2013, 2014) and through a dense network of static cells (Lee et al. 2017).In particular, with respect to previous mechanical models (Giverso et al. 2014(Giverso et al. , 2018)), this work provides better insights into the phenomenon, since the whole dynamics of the process is investigated, along with the influence of the cell membrane.Furthermore, the motion of the cell in this case does not require the presence of an external flux as in Lee et al. (2017), but only relies on cell deformation, cortex polymerization and the microtubule activity.Finally, compared with previous models derived using an extended version of a Cellular Potts Model (Scianna et al. 2013;Scianna andPreziosi 2013, 2014), this work allows, in principle, to obtain quantitative results of the whole migratory process, by including into the model identifiable mechanical parameters.
Of course, the dynamics of cell and nucleus deformation and motion that can be reproduced by our model is wider than the one captured by these benchmark simulations, therefore we will investigate in Sect.4.3 the dependence of the mean cell speed and the nucleus shape on the parameters of the model.

Influence of the Model Parameters
In order to understand how the different parameters of the model affect the capability of the cell to migrate inside sinusoidal channels and its mean speed, we perform a set of numerical experiments, changing one parameter at a time.Indeed, the model put forward in the present paper makes it hard to predict analytically the average velocity of a cell moving inside a structured channels.Therefore, we should rely on numerical simulations.To compare the results, we defined the mean cell speed as the average speed of the tip of the cell over one period of the cyclic motion, and the mean nucleus area as the average area enclosed by the nuclear membrane over the same period.The value of the mean cell speed for varying nuclear stiffness is shown in Fig. 6a, with the area-change constraint (μ n ) and the bending modulus of the nucleus (k b ) along the horizontal and vertical axes, respectively.From Fig. 6a, one sees that the capability of the cell to move inside the channel and its speed highly depends on the nucleus mechanical properties that determine nucleus ability to deform, as it is apparent from the insets.These show the deformation of the nucleus at a given instant of time for different sets of parameters (marked with coloured dots in the parameter space).
Indeed, cell migration is hampered when the energy needed to bend the nuclear envelope is too high (i.e. for high values of the parameter k b , corresponding to the red dot in the parameter space), since in this case the nucleus is not able to deform much and acquire the hourglass shape required to pass the constrictions in the sinusoidal channel and it stays round, occupying the space between two constrictions.Conversely, decreasing the value of k b (see the gray dot in Fig. 6a and the inset with the corresponding nucleus shape), and keeping the value of μ n fixed, the nucleus can deform and pass through constrictions.Therefore, for a given value of μ n , the speed of the cell decreases for increasing values of the bending modulus k b , until the cell is stuck inside the channel (zero speed).The threshold for k b allowing cell motion depends on the value of μ n .
Specifically, the constraint imposed by the bending stiffness of the nucleus is more restrictive when the chromatin inside the nucleus is little compactible, i.e. for high values of the parameter μ n .Indeed, when μ n is sufficiently high, the nuclear deformations occur while maintaining the nuclear area close to the target one A * n .This is evident looking at the plot in Fig. 6b, where we report the average nucleus area (still computed over one period of the cyclic motion) with respect to the nucleus target area, for the same values of the nucleus mechanical parameters used in Fig. 6a: for high values of the parameter μ n the nucleus moves preserving in average more than the 90% of its target area.In this case, the deformed nucleus even pinches twice (see the deformed nucleus in correspondence of the yellow dot in Fig. 6).
The plot in Fig. 6b also allows to comment on the admissibility of the cell velocities predicted by the model.Indeed, from Fig. 6b it is clear that for small values of the parameters μ n , the nucleus can decrease its area (which would correspond to the volume in a three-dimensional simulation) below a physiological threshold (shaded region on the left of Fig. 6b).Even though there is biological evidence (Friedl et al. 2011;Rowat et al. 2006;Versaevel et al. 2012) that the volume of the nucleus is not preserved during large elongations-which suggests that the nuclear envelope is permeable to aqueous material and that the chromatin structure can compact itself (chromatin condensation)-the nucleus volume cannot shrink under a minimum threshold.In particular, in Rowat et al. (2006), it was shown that although nuclei experienced a marked loss of total volume under aspiration, it stabilized above 30-40% of the initial nuclear volume.Therefore, mechanical parameters allowing the nuclear area going below the 30-40% of the target area should be disregarded (see the white isoline in Fig. 6 Average cell speed (a) and ratio between the average and the target nucleus areas (b) while moving inside a sinusoidal channel.The shaded region on the left corresponds to a ratio A/A * n < 50%.The region bounded by the dotted line corresponds to A/A * n > 70% and zero speed.The black dot denotes the point of maximum speed, while region for which the velocity is above 99% of the maximum velocity is marked with the dashed white line.The blue dot corresponds to the values of μ n and k b taken for Fig. 7.For select parameters, marked by coloured dots, the cell and the nucleus are illustrated.The parameters used are given in Table 1 (Color figure online) Fig. 6a, corresponding to the threshold A/A * n = 50%).We observe that, relaxing the constraint on the area change (i.e.lowering the value of μ n ), so that the nucleus can shrink, the cell can move inside the channel even for high values of bending modulus k b , albeit with lower velocity (green region in the top of Fig. 6a).In this case, the nucleus maintains an elongated ellipsoidal shape and does not acquire an hourglass deformation as the cell moves inside the channel (see the nuclear shape corresponding to the white dot in Fig. 6).
Furthermore, the proposed model predicts the existence of an optimal region in the space of the mechanical parameters k b -μ n for which the cell speed is maximal (region delimited by the white dashed line in Fig. 6a, with the black dot corresponding to the maximum cell speed) and the nuclear area is in the range 70-85% of the target area.
Finally, in Fig. 7, we consider the cell average velocity for varying values of the parameters describing the channel geometry and we compare the results obtained either with or without the description of the nucleus.In particular, we first consider the pulsation of the channel in determining the cell capability to move.As observed for the cell without the nucleus (blue line in Fig. 7a) the cell is able to migrate only Nucleus bending stiffness Nucleus/cortex interaction charact.length ξ cont 10 Nucleus/cortex interaction coefficient k cont 5 Centrosome related parameters Centrosome link stiffness k e 10 −3 if the channel is sufficiently structured, i.e. for channel pulsation above a minimum threshold, f ω min .

Channel geometry related parameters
The value of the threshold f ω min is only slightly influenced by the presence of the nucleus.The latter mainly affects the cell average speed, which decreases in the presence of the nucleus.Specifically, the cell speed decreases when increasing the value of k b (see the yellow and red curves).However, for high values of k b , the cell gets stuck (see the red curve in Fig. 7a) even for intermediate values of channel pulsation, since the hourglass nuclear deformation is impeded by the high bending modulus and the space between two subsequent constriction is enough to host the cell nucleus.Increasing the channel pulsation further, the nucleus characterised by a high  1 bending modulus is again able to migrate since the space between two constriction becomes restrictive for the preservation of the nuclear area and the nucleus can deform acquiring an ellipsoidal shape, without following the structure of the channel walls.We remark that, since the discretisation of the cell membrane is finite, the channel structure will not be well resolved as f ω 0 is taken larger and larger.This explains why the first plot in Fig. 7a is cut on the right hand side.Note that this will occur for any discretisation size.
The influence of the channel width on the cell velocity is analysed in Fig. 7b: without considering the nucleus (blue curve), the cell can move even for smaller channel width and this "physical limit of cell migration" is related to the capability of the cell envelope and the cytosol to deform and enter even small gaps.Indeed, the size of the smallest constriction inside the channel is 2( f width − f β ).When we include the nucleus, the cell is no longer able to move inside channels with a neck which is too narrow (i.e.small values of f width ), since the nucleus hinders cell capability to deform (yellow curve).Raising the bending modulus (red curve), the cell can move inside the channel for large values of f width only.On the other hand, the upper limit is the same in all three cases, since it is related to the capability of the whole cell to maintain the contact with the channel wall and it is not affected by the presence and the mechanical properties of the nucleus.Furthermore, the plots in Fig. 7b show the well known bimodal behaviour of cell velocity for varying channel width (Ulrich et al. 2009;DiMilla et al. 1993;Kuntz and Saltzman 1997): cells cannot migrate both inside very small channels (since deformations would be too large) and inside very large channels (since the cell does not touch the boundary of the channel.The velocity of cell motion is also in this case slowed down by the presence of the nucleus.Specifically, the speed decreases for increasing values of k b , as previously observed. The bimodal behaviour in the cell velocity can be observed also when changing the value of f β and keeping the channel width fixed (see Fig. 7c).In this case, the upper limit of f β is related to size of the small constriction that the cell could enter in order to move inside the channel and thus it is highly influenced by the presence of the nucleus and by its mechanical properties.On the other hand, when f β is small the channel is almost flat and the cell cannot migrate, independently on the presence of the nucleus and its mechanical properties.Therefore, the lower limit for f β is the same both for the cell with and without the nucleus, since it is related to the mechanism of motion, which requires a sufficient structure of the lateral walls.

Conclusions
In this work we have proposed a continuous mechanical model describing single cell adhesion-independent migration inside restrictive 3D environments.The cell migration is driven by a local imbalance in the polymerization and depolymerization of the actin network underneath the cell membrane, which induces a cortex flow.The cell shape is determined by the balance of cytoplasmic pressure, elastic behaviour of the cell cortex and interactions with the subcellular elements and the channel walls.Differently from previous works (Jankowiak et al. 2020), the influence of the nucleus in the process of cell migration is explicitly taken into account by including an inner nuclear membrane, connected to the cell membrane by the microtubule structure, responsible for the positioning of the nucleus inside the cell.The nucleus shows a resistance to bending, to stretching and to changes in the area enclosed by the nuclear membrane and it becomes the limiting factor in determining the ability of the cell to migrate across small neckings of the channel.Therefore, the proposed model represents an advancement with respect to the state of the art, since it provides a purely mechanical description of adhesion-free migration, without requiring external chemical stimuli as done in Lee et al. (2017), Chen et al. (2018) and Cao et al. (2016).It also allows testing the influence of nucleus mechanical properties in the determination of the physical limit of cell migration, which has been neglected in previous models (Wu et al. 2018;Stotsky and Othmer 2022;Kaoui et al. 2008;Moure andGomez 2017, 2018;Jankowiak et al. 2020).
The model equations have been discretised and solved numerically in order to simulate the process in a 2D geometry, corresponding to a section of the 3D channel with lateral structured walls and top and bottom flat walls.To formulate an efficient and stable (under appropriate time step restrictions) discretisation scheme, we adapt the approach proposed in Mikula and Ševčovič (2004) and Beneš et al. (2009) to our model.
The numerical simulations reproduce qualitatively the behaviour observed in the biological experiments (Reversat et al. 2020), where it is pointed out that adhesionindependent migration needs both confinement and sufficiently structured channel walls.This behaviour is purely related to the cell membrane behaviour and it is not influenced by the presence of the nucleus, since the threshold for the channel structure depth (i.e.lowest value of f β that allows cell motion) predicted by our model is the same of the one predicted by the model developed in Jankowiak et al. (2020).However, the cell speed inside the channel and the physical limit of cell migration (i.e. the size of the smallest opening in the channel) substantially depends on the presence of the nucleus and on its mechanical properties, as demonstrated by the parametric study showing the dependency on geometric properties of the channel and on the mechanical properties of the nucleus.In particular, in agreement with biological experiments (Wolf et al. 2003(Wolf et al. , 2013)), we observe that a little-deformable nucleus is a limiting factor for cell migration inside restrictive environments.Indeed, by keeping both the geometry of the channel and the cell membrane parameters fixed, we can show the transition from a migrating cell to a non-moving cell by changing the mechanical properties of the nucleus only (i.e.increasing the bending modulus k b and the elastic area constraint μ n ).
Furthermore, the sensitivity analysis of the mechanical parameters of the nucleus predicts the optimal nuclear deformability, for which the cell reaches its maximal speed, for low values of k b and intermediate values of μ n .
Without claiming to provide quantitative numerical measurements, the results presented in this work are a proof of concept for a comprehensive model of single cell adhesion-independent migration taking into account cell and nucleus mechanics.Nonetheless, these results need to be validated quantitatively by comparing the predicted evolution with the actual spatio-temporal evolution of the moving cell inside a structured channel, by also keeping track of nuclear deformations.Therefore, additional work is needed to investigate more realistic in vitro and in vivo conditions, in order to quantitatively validate the model.
From a modelling point of view, this study has to be seen as a first step and several components would benefit from a more detailed treatment.In particular, the role of the compensating force, that can be attributed to the transport of actin monomers and their polymerization, certainly requires further studies.For example, one could include an explicit description of the acto-myosin machinery and its influence in cell polymerization.Another interesting direction would be to take into account the resistance to deformation of the cortex, which can then be considered as visco-elastic and not only elastic.Moreover, the mechanical description of the nucleus is kept rather simple in the present model and a detailed description of both the deformations occurring inside the nucleus and the exchange of liquid between the nucleus and the surrounding cytosol during the whole process of cell migration is still missing.Finally, the microtubule description is kept really simple and it does not take into account some important biological observations, such as the bending of some MTs close to the cell cortex, the dynamic instability process, the detailed description of the anchorage of the MT structure to the cell cortex and to the intermediate filaments network surrounding the nucleus membrane.The inclusion of these biological observation inside the model is fundamental to give a biologically sound description of the active and passive force exerted by the MTs and of their role in positioning the cell nucleus, while maintaining the cell shape.From the biological point of view, it would be interesting to perform ad-hoc experiments to quantitatively verify the model prediction.Indeed, even though most of the model parameters could be in principle measured or at least estimated from numerical experiments, it is not possible to derive all of them from the biological observations reported in the literature.Therefore, future works will certainly be addressed to the fitting of the model parameters with experimental data.
Taking all these effects into account could lead to a more comprehensive understanding of the multiple factors involved in determining the physical limit of cell migration during non-proteolytic migration of cells.

Appendix A: Nuclear Energy Variation
We here compute the variation δ E n of the nuclear membrane energy function, in order to obtain the force acting on the nuclear membrane.We named In order to compute δ E ( j) n , with j = 1, 2, 3, 4, it is convenient to use the parametrization of the curve with respect to σ , since d will undergo variations as well for displacement of Y .
The first energy term with respect to σ (assuming K 0 = 0) reads and its variation is (Kaoui et al. 2008) The variation of the second energy term can be easily computed In order to calculate δ E (3) n and δ E (4) n , we use the identity Finally, the nuclear membrane mechanical force reads where X c is the centrosome.Then the line segment between X i and X i+1 has the following equation in polar coordinates where Then we have the following identities: In order to integrate function over the cortex f ( MT (θ )) = f (s MT (θ )) w.r.t.θ , we need to know MT (θ ) and −1 MT (s) in order to perform the change of variables The expression of | −1 MT (s)| is derived from (7): Therefore, we can define the local length element g:=|∂ σ Y (σ )|, i.e. d = gdσ The curvature is denoted by K and ν is the tangential angle, so that T = (cos ν, sin ν) and N = (sin ν, − cos ν).The region enclosed by n is called n and d A denotes its surface measure.The curve's total energy E n is given by: where the first three terms are given by (A1), while W denotes a potential, which combines the contribution related to the nuclear membrane surface tension (i.e.E (4) n given by (A1)), the contact interaction with the cell cortex and the elastic constraint linking the nucleus with the centrosome.Following the method proposed by Mikula and Ševčovič (2004), the evolution of the planar curve n is, then, given by: where K , ν and g = |∂ σ Y | are the curvature, tangential angle and local length element, respectively.Keeping in mind the previous governing equations, it is possible to calculate the time derivative of the energy functional E, formally given by Then, the choice of β is made to maximise the decrease of energy due to the normal component, whereas the tangential velocity α, in principle, is a free parameter.In this paper, we refer to the asymptotically uniform tangential redistribution derived in Mikula and Ševčovič (2004) and Beneš et al. (2009), to impose the evolution of α in order to avoid nodes concentrating on points, which would lead to poor approximation and eventually inversion of ill-conditioned matrices.Therefore, we take where .= L −1 n n .d denotes the average over the whole curve n , which makes (36) a non-local equation, whereas ζ > 0 is a given positive constant.
Recalling the identity: so that, setting g = exp η, we can rewrite (35) as We remark that the second equation in (38) has been obtained remembering that K = ∂ ν.

C.2 Quadrature Formulae for the Mikula Formulation
We can now write the corresponding discretised equations of (38).To do so, we discretise the fixed parameterization interval [0, 1] uniformly in N 2 subintervals, each of equal length h = 1/N 2 and indexed by i ∈ {0, N 2 − 1}.Time is also discretised with a time-step t.The discrete time is indexed by j ∈ N, so that t j = j t.In what follows, indices corresponding to space (resp.time) will be in subscript (resp.superscript), so that the point Y (ih, j t) is written

2
. The length of the dual finite element is denoted by q j i = 1 2 r j i + r j i+1 .Then, the system of equations ( 36)-( 37)-( 38) is solved for the discrete quantities α In particular α and we discretise in time to get where L j = l r j l , B j = (L j ) −1 l r j l K j l β j l and β j i obtained by the discretization of (37), i.e.
This determines α j i for all j, provided a closure for α j 0 .It is possible to assume α j 0 = 0, i.e. the point Y j 0 moves along the normal (Beneš et al. 2009).The local length of the finite element is updated in a similar way, using (38) c : If we approximate the time derivative as t −1 (K ), K j i solves a linear system of the form (40) Then, we have to solve the tangent angle equation (38) b , using the following approximation and approximating the time differential by finite differences, we obtain the following system for ν So that the complete semi-discrete equation is: Again approximating the time differential by finite differences, we get a system of the following form for Y j i : The coefficients of the linear equations ( 40), ( 42) and ( 44) are reported in the following.

C.2.1 Discretisation of the Equation for K
In the following, we report how the different terms appearing in Eq. ( 39) have been discretised Therefore, the coefficients of ( 40) are given as follows:

C.2.2 Discretization of the Equation for
In the following, we report how the different terms appearing in Eq. ( 41) have been discretized Then, the coefficients of (42) are ,

C.2.3 Discretisation of the Equation for Y i
In the following, we report how the different terms appearing in Eq. ( 43) have been discretised The coefficients of (44) are We remark that the linear dependency on Y inside W (Y ) can be possibly treated implicitly.

Fig. 2
Fig.2Left: Parameterization of the cell at time t.The front and back of the cell are shown, with the corresponding polymerization regions for a possible choice of f .The dark yellow and light yellow area highlight the support of f + and f − , respectively.Right: Parameterization at time t + ε.The dashed lines show how the material points at t and s = 0 and s = s 1 are mapped at t + ε.Because of the polymerization, ∂s/∂t > 0 in the upper part of the curve between the two yellow regions, so that, at t + ε, the material point which was located at X (t, s 1 (t)) is at X (t + ε, s 2 (t + ε)), for some s 2 > s 1 (Color figure online) b a φ ( f (x)) dx for any convex function φ : R → R and f ∈ L 1 (a, b), see e.g.Evans (2010, B, Theorem 2).

Fig. 3
Fig.3Left: The sinusoidal channels used in the numerical simulation are described by the parameters: (1) f width , representing half of the mean width of the channel, (2) f β , which is the amplitude of the oscillation of the wall, (3) f ω 0 , which is the pulsation of the sinusoidal channel.Right: Illustration of the construction of the initial data C 0 and its components C 0 i .The value of x min 0 is fixed to −π/2 f ω 0 , and x max 0 is chosen so that the initial cell area (blue) is equal to A * c (Color figure online)

Fig. 4
Fig. 4 Equilibrium configuration for f β = 0, the values for the remaining parameters are presented in Table 1.The cell is in blue and the nucleus in light blue.The dark blue dot is the centrosome.The thickness of the yellow (resp.light yellow) region on the cell membrane indicate the strength of the polymerization (resp.depolymerization).The arrows indicate the flow of the cortex relative to the cell (Color figure online) b = 10 −2.5 , μ n = 50); (c) the migrating cell with an intermediate bending modulus of the nucleus (Fig. 5c, with k b = 10 −1.5 , μ n = 50); (c') the migrating cell with a higher value of the relaxation parameter μ n , and the same bending modulus as in (c) (Fig. 5c', with k b = 10 −1.5 , μ n = 100); (d) the non-migrating cell with a high bending modulus of the nucleus (Fig. 5d, with k b = 10 −0.5 , μ n = 50).

Fig. 5
Fig. 5 Snapshots of the cortex/nucleus/centrosome system at different times, without nucleus (a) and with increasing nucleus stiffness k b .The average velocity decreases in the presence of the nucleus, and decreases further as k b increases [k b = 10 −2.5 (b), k b = 10 −1.5 (c)], eventually reaching 0 [k b = 10 −0.5 (d)].The row (c') illustrates the situation for a larger value of the nucleus area constraint relaxation parameter μ n = 100.Other rows correspond to μ n = 50.The remaining parameters are presented in Table 1

Fig. 7 (
Fig. 7 (Color figure online) Comparison of the motilities associated with cell with and without nucleus, for three different channel parameters.Top: Wave number of the channel pattern f ω 0 .Middle: mean width of the channel f width .Bottom: depth of the pattern f β .Model parameters are summed up in Table 1 The dependent variables W (Y ) and ∇W (Y ) are consequently evaluated numerically at nodes, i.e.W j i = W (Y j i ) = W (Y (ih, j t)), whereas the tangential and normal vectors are approximated on the finite element [.We first integrate (36) on the finite element [Y i−1 , Y i ] discretisation of the equation for Y (38) d is obtained by integrating on the dual element [ Ỹi+1 , Ỹi ], i.e.

Table 1
Values and ranges of the dimensionless parameters used for the numerical experiments Integrating Eq. (38) a on [Y i−1 , Y i ] we get j i ).