Nuclear Stress-Strain State over Micropillars: A Mechanical In silico Study

Cells adapt to their environment and stimuli of different origin. During confined migration through sub-cellular and sub-nuclear pores, they can undergo large strains and the nucleus, the most voluminous and the stiffest organelle, plays a critical role. Recently, patterned microfluidic devices have been employed to analyze the cell mechanical behavior and the nucleus self-deformations. In this paper, we present an in silico model to simulate the interactions between the cell and the underneath microstructured substrate under the effect of the sole gravity. The model lays on mechanical features only and it has the potential to assess the contribution of the nuclear mechanics on the cell global behavior. The cell is constituted by the membrane, the cytosol, the lamina, and the nucleoplasm. Each organelle is described through a constitutive law defined by specific mechanical parameters, and it is composed of a fluid and a solid phase leading to a viscoelastic behavior. Our main objective is to evaluate the influence of such mechanical components on the nucleus behavior. We have quantified the stress and strain distributions in the nucleus, which could be responsible of specific phenomena such as the lamina rupture or the expression of stretch-sensitive proteins.


Introduction
Cells continuously adapt themselves to their environment and stimuli they received (i.e., chemical, electrical, mechanical, …) [1][2][3]. More specifically, they are able to undergo large strains during confined migration through sub-cellular or sub-nuclear pores. During such a process, the nucleus, the most voluminous and the stiffest cellular organelle, plays a critical role [4,5]. Some cells such as cancerous cells can even undergo the rupture of the nuclear lamina to be able to migrate across healthy tissues [6,7]. Therefore, quantifying nucleus strains and stresses can be crucial to diagnose cancer and other pathologies in patients.
To do so, patterned microfluidic devices have been employed during the last few years in order to characterize the cell mechanical behaviour [8,9] and the nucleus self-deformations [10][11][12][13] and shape changes [14,15] induced by mechanical forces, which are due to the interaction between the cell and the topological surface. Assays on micropillared substrates involve successive steps: (i) contact between the cell and the pillars, (ii) adhesion of the cell on the pillars surface, (iii) cell spreading, (iv) cell polarization and (v) cell crawling.
A series of analytical and numerical models exist in the literature focusing on the interactions between the cell and a flat substrate. The former provides information on the spreading process with no excessive computational cost [16][17][18][19]. The latter, which may be discrete [20][21][22] or continuum [15,[23][24][25], allows to investigate the intra-cellular rearrangement or to obtain quantitative results at the global or the local scale. In our previous paper [26], we have proposed an in silico two-dimensional (2D) model which simulates the first three steps (i.e., contact, adhesion and spreading) of the interaction between the cell and the micropillared substrate and provides insights on the mechanisms inducing nuclear deformation. We have been able to determine the role of the gravity and of the actin fibers above and beneath the nucleus responsible for a pushing and a pulling force, respectively.
In the present paper, we have adapted the model presented in [26] and we have focused on step one only (i.e., contact between the cell and the pillars). The cell and nucleus behaviours are described using specific mechanical tools (i.e., constitutive laws, mechanical properties, fluid and solid phase balance). By performing a sensibility study, our objective has been to determine the influence of these parameters on the interaction between the cell and underneath micropillared substrate. Then, the model provides the nucleus stress and strain fields over time and gives insights on the global cellular behavior that can be further explored experimentally.
In Section 2, we describe the mathematical framework of the model including the geometry (Section 2.1), the constitutive laws governing the behaviour of the cell and of its components (Section 2.2), the external forces applied to the cell (Sections 2.3 and 2.4) and the numerical implementation (Section 2.5). The results of the different simulations are presented in Section 3 and some conclusions and perspectives are drawn in Section 4.

Cell Geometry
Given the symmetry conditions, we consider here the cell c in its initial configuration as a semicircle of radius r c (Fig. 1). It is constituted by the membrane m (external radius r m ¼ r c and thickness t m ), the cytosol cs (external radius r cs ), the lamina l (external radius r l and thickness t l ) and the nucleoplasm np (external radius r np ) (Fig. 1). The membrane and the cytosol form the cytoplasm cp , whereas the lamina and the nucleoplasm form the nucleus n . Each component of the cell i is described via a spatial characteristic function g i , which is the composition of a regularized Heaviside function H (equal to 0 and 1 for negative and positive argument, respectively), and a level set function l i . Then, the characteristic functions describing the cell and its components are the following: g cs ¼ H l cs À g l À g np ¼ H kp À c 2 c;p k À r 2 cs À g l À g np (4) g m ¼ H l m À g l À g np À g cs ¼ H kp À c 2 c;p k À r 2 m À g l À g np À g cs (5) with k Á k the euclidean norm of a vector, p the initial position of any point of the system and c c;p the center of the cell of coordinates c x ; c y À Á .

Constitutive Laws
Each cell component is composed of a solid and a fluid phase. To simplify the approach, we assume that both phases deform in parallel like in a Kelvin-Voigt model. Consequently, the overall stress S and the deformation F can be expressed as where the subscripts s and f indicate the solid and the fluid phases, respectively, while the subscript i indicates a specific constitutive law for the solid phase, as described in the followings (Section 2.2.1), and c s and c f the respective concentrations are given by c s ¼ c s;m g m þ c s;cs g cs þ c s;l g l þ c s;np g np (8)

The Solid Phase
When describing the solid phase of an isotropic elastic material, for stability reasons, its strain energy must be poly-convex with respect to the three invariants I 1 , I 2 and I 3 [27,28], which are expressed as follows: From a physical point of view, I 1 describes the length variations of the system in the three directions, I 2 gives both the shear and length deformations and I 3 describes the volume variations of the system [29,30]. Then, given the multistructural organization of the cell, it is important that all the invariants are considered. Nonetheless, in the present model, four different materials have been tested to point out their specific features. Specifically, we have implemented: 1. A standard Saint-Venant material, which only depends on the first and second invariants I 1 and I 2 [31,32] and it is able to capture the deformations along the lines (i.e., the fibers deformations); 2. A Neo-Hookean compressible material, which depends on the first I 1 and third I 3 invariants [33,34] and allows to take into account the deformations along the lines and inside the volumes (i.e., both fibers and cytosol plus nucleoplasm deformations); 3. The Mooney-Rivlin and Yeoh compressible materials which depend on the three invariants [25,35] and describe the deformations along lines, volumes and surfaces (i.e., the membrane and the lamina). Furthermore, the Yeoh material takes into account the successive softening and stiffening behavior of an elastic material made of fibers.
The second Piola Kirchhoff stress S s;i , where the subscript i indicates the material, is expressed as follows: for the standard Saint-Venant material S s;SV reads for the compressible neo-Hookean material S s;NH is given by for the Mooney-Rivlin material S s;MR can be written as -finally, for the Yeoh material S s;Y reads where a Y 1 ¼ À0:5a MR and a Y 2 ¼ 0:1a MR are two constants [26].

The Fluid Phase
For the fluid phase of the cell, a Newtonian viscous fluid is considered but it must be defined in the Lagrangian configuration in order to ensure the compatibility with the solid phase. Thus, the Cauchy stress r f is classically given by where k f and l f are the isotropic and deviatoric viscosities, respectively and D f must be expressed as a function of the rate of the dilatation tensor C f as follows: with the superscript T indicating the transpose of a matrix. (7)), substituting r f in Eq. (21), we obtain the expression of the second Piola-Kirchhoff

The Gravity Force
The cell is submitted to the gravity force f g f g ¼ Àq c t g gi y (24) where q c is the initial cell density, g the gravitational acceleration and i y the vertical unit vector. t g is employed here to smoothly apply the gravity force to ensure the computation fast convergence. It is the composition of a regularized Heaviside function H and a level set function l g and it reads with t the time and T g a constant.

The Cell Environment
As the gravity is applied, the cell settles down and starts interacting with the underneath micropillared substrate, which is constituted by a micropillar and a flat region and it is defined by a characteristic function g s as follows: where l mp and l flat are two level set functions, x and y are the actual coordinates of any point in the system, x mp and y mp define the position of the micropillars, w mp and h mp define the width and the height of the micropillar, respectively and y flat is the position of the flat region with respect to the y-axis.
Once the cell approaches the substrate, the contact force f ct applies over a very thin layer corresponding to the intersection between the cell and the substrate. We have shown that f ct can be approximated by a volume force [26] and it can be computed via a penalization technique as follows: where l ct is the penalization coefficient. cof F ð Þ is the cofactor matrix defined as J c F ÀT and it is employed here to bring back to the initial configuration the outward normal vectors to the micropillar (n mp ) and to the flat substrate (n flat ), which are obtained as follows:

Finite Element Model
The global equilibrium of the system in the initial configuration can be expressed as with Div the first order divergence operator in the initial configuration, and a the acceleration.
In order to employ a classical finite element approach, we multiply each term of Eq. (30) by the kinematically admissible displacement test function w and we integrate over the cellular domain c . Then, by using some algebraic operations and by applying the Stokes theorem, we obtain where n c is the outward normal to the cell and the first and the third term indicate the internal cellular stress and the volume forces, respectively. The second term describes the surface forces applied to the cell and specifically the contact force f ct between the cell and underneath micropillared substrate. As mentioned above (Section 2.4), f ct is applied through a very thin layer. Thus, the surface integral in Eq. (31) can be written as a volume integral over the penalization thickness h p . Therefore, Eq. (31) becomes Eq. (26) has been manually implemented using the weak form tool in COMSOL Multiphysics. The spatial discretization is obtained via quadratic polynomials for each isoparametric element of the mesh (mesh size between 0.3 µm and 1 µm). The time discretization is achieved via a second-order backward differentiation formula (BDF). The solution is computed using a nonlinear Newton scheme with a relative tolerance of 1% on the displacement error estimation.

Results and Discussion
The radii of the cell (r c ) and the nucleus (r n ) have been set equal to 10 µm and 5 µm, respectively. The membrane and the lamina have a thickness t m and t l of 0.5 µm and 0.1 µm [36], respectively (Fig. 1). The cell density q c and the Poisson ratio t c have been set to 1000 kg/m 3 and 0.45, respectively. The isotropic (k f ) and deviatoric (l f ) viscosities are equal to 1000 Pa·s and 5 × 10 −3 Pa·s, respectively. The Young's moduli have been set to E m = 100 Pa [31], E cs = 10 Pa [37], E l = 1000 Pa [31,38] and E np = 10 Pa [39]. The gravity force is gently applied on the cell and reaches its maximum at 2T g ¼ 10000 s. Through the gravity, the cell comes in contact with the micropillar which has a width (w mp ) and a height (h mp ) of 20 µm and 9.5 µm, respectively (Fig. 1). All the simulations last for T end = 18000 s. All the parameters of the model are listed in Table 1. For each series of simulations that we are going to present, the nucleus stress has been evaluated in terms of Mises Cauchy stress r Mises , whereas the strains have been quantified through three variables: (i) J 1 =I 1 = Tr C ð Þ to detect the fibers elongation, (ii) J 2 = Tr C 2 À Á to evaluate the shear and length deformations and (iii) J 3 = I 3 = det C ð Þ to assess the surface variations. End of the simulation 18000 s

Cell Components Constitutive Law
For the first series of simulations, we have tested different constitutive laws for the cell components. All the parameters defining each constitutive law are reported in Table 1. The fluid coefficients have been set equal to c f ;m ¼ c f ;cs ¼ c f ;l ¼ c f ;np ¼ 0:5.
The Saint-Venant material model is the simplest and mostly used in the literature, but it is not robust enough to describe very large deformations. The Neo-Hookean and the Mooney-Rivlin materials are efficient in considering large strains, but they do not exhibit typical stiffness-softening followed by hardening during deformation. Finally, the Yeoh model is theoretically the most consistent, but it is much more complex than the others due to a higher number of parameters to define (Eq. (21)).
The evolution as a function of time of the Mises stress r Mises and J 1 , J 2 and J 3 in the nucleus n are shown in Fig. 2. The highest value of the von Mises stress (1.042 Pa) is obtained for the Mooney-Rivlin material (Fig. 2a). The Neo-Hookean and the Yeoh materials are very close with a maximum of 1.032 Pa and 1.027 Pa, respectively. For the Saint-Venant material, the maximum is slightly lower (0.86 Pa). We can estimate the contribution of the nucleus to the total cell deformation. In terms of fibers elongation (J 1 ), we found that when implementing the Mooney Rivlin material the nucleus undergoes 56% of the global deformation, whereas with the Saint-Venant material only 49%. Regarding the shear and length variations (J 2 ), the Mooney Rivlin ensures again the highest nucleus contribution (81%) vs. the Yeoh material for which the nucleus supplies for only 61%. Finally, for the surface variations (J 3 ), contributions between 18% (Saint-Venant material) and 23% (Mooney Rivlin material) are found. In Fig. 3 we show the overall deformation of the cell and the nucleus at t = 18000 s for the different material models.
Since the cell is composed by several organelles having different structures and geometries, there is a need for a constitutive macroscopic law able to describe the stresses associated to such heterogeneous strains. More specifically, taking into account these stresses implies low values of J 1 , J 2 and J 3 . Thus, according to our results, the Yeoh material model seems to be the most appropriate, which confirms the theoretical remarks at the beginning of the section.

Nucleus Fluid and Solid Phases
For the second series of simulations we have evaluated the influence of the fluid-solid phases of the lamina and of the nucleoplasm on the nucleus behaviour. Specifically, we have let vary c f ;l and c f ;np from 0.1 to 0.9, testing all the combinations for a total of 24 simulations (one combination has not converged, i.e., c f ;l ¼ c f ;np ¼ 0:9). The values of c f ;m and c f ;cs have been set to 0.5. The Young moduli of the cell organelles are those defined at the beginning of Section 3.
In Fig. 4, values of r Mises (Fig. 4a) and of J 1 ; J 2 and J 3 (Figs. 4b, 4c and 4d) in the nucleus n are presented. r Mises is maximal for c f ;l ¼ 0:1 and c f ;np ¼ 0:7 and it is equal to 1.32 Pa.
On the one side, the maximal strains are found for c f ;l ¼ 0:7 and c f ;np ¼ 0:9 and they are equal to J 1 ¼ 4:08 (71% of the total cell deformation), J 2 ¼ 15:20 (163% of the total cell deformation) and J 3 ¼ 1:31 (20% of the total cell deformation). This combination corresponds to r Mises ¼ 0:45 Pa. In Fig. 5a, the total deformation of the cell organelles at t = 18000 s is shown and one can notice that the nucleus is particularly stretched. Such a large strain could exacerbate the sensor functions of the nucleus as it has been observed in previous experimental studies [40][41][42] and therefore induce cell motility via the activation of specific stretch-sensitive proteins.
On the other side, the lowest values of J 1 and J 2 (4:08 and 15:20, respectively) are found for c f ;l ¼ 0:3 and c f ;np ¼ 0:1, with r Mises ¼ 1:08 Pa, and they correspond to 51% and 58% of the total cell deformation. Whereas the minimal value of J 3 (1.31) has been obtained for c f ;l ¼ 0:9 and c f ;np ¼ 0:1, with r Mises ¼ 0:85 Pa, and it corresponds to 23% of the total cell deformation.

Nucleus Mechanical Properties
For the last series of simulations, we have let vary the Young's moduli of the nucleoplasm (E np ) and of the lamina (E l ). Specifically, we have set E np = 10, 30, 50, 70, 90 Pa and E l = 1000, 1500, 2000, 2500, 3000 Pa for a total of 25 simulations, which have all converged. For the fluid and solid phases of the organelles, we have set here c f ;np ¼ c f ;l ¼ c f ;cs ¼ c f ;m ¼ 0:5. Fig. 6, we show the results in terms of Mises Cauchy stress r Mises (Fig. 6a) and of J 1 ; J 2 and J 3 (Figs. 6b,  6c and 6d) in the nucleus n . r Mises is maximal for E np = 10 Pa and E l = 2000 Pa and it is equal to 1.47 Pa. The maximal values of J 1 and J 2 are observed for E np = 10 Pa and E l = 1000 Pa (r Mises ¼ 1:03 Pa) and they are equal to 2.44 (53% of the total cell deformation) and 4.27 (69% of the total cell deformation), respectively) (Fig. 5b). For J 3 , the maximal value is equal to 1.16 (22% of the total cell deformation) and  it is found for E np = 10 Pa and E l = 3000 Pa, with r Mises ¼ 1:47 Pa. As for the minimal values of J 1 and J 2 , they are equal to 2.11 (50% of the total cell deformation) and 2.26 (51% of the total cell deformation), respectively. They are found for E np = 90 Pa and E l = 1500 Pa and provide r Mises ¼ 0:55 Pa. The lowest J 3 has been found equal to 1.10 (24% of the total cell deformation) for E np = 90 Pa and E l = 1000 Pa with r Mises ¼ 0:46 Pa.

Comparison with Experimental Data
In the previous series of simulations, t end has been set equal to 18000 s. In this section we propose to go further (t end = 24000 s) in order to explore more severe nucleus deformations and to be able to compare our numerical results with experimental observations. We have run three additional simulations (S 1 , S 2 , S 3 ) for which a Saint-Venant material model has been implemented. However, the mechanical properties of the nucleus (i.e., E np and E l ) as well as the nucleus fluid concentration (i.e., c f ;np and c f ;l ) have been modified as reported in Table 2.
From a quantitative point of view, the results in terms of maximal stress (r Mises ) and strain (J 1 , J 2 , J 3 ) in the nucleus n are reported in Table 3.
In Fig. 7, the total cell and nucleus deformation at t = 24000 s is shown. One can notice that for S 1 (Fig. 7a), the final shape of the nucleus is not significantly different from the previous simulations. However, as the fluid concentration of the nucleoplasm and the lamina is increased (S 2 ), the nucleus undergoes a higher deformation (Fig. 7b) leading to an important reduction of the distance between the nuclear lamina and the micropillar boundary. Additionally, the nucleus starts acquiring a 'peanut' shape  as it has been experimentally observed for different cellular phenotypes including mesenchymal stem cells, osteosarcoma cells and bone marrow stroma cells [11][12][13][14]. Such an outcome becomes even more evident when the nuclear lamina is ablated (S 3 ) (Fig. 7c), which correspond to a decrease of the lamina Young's modulus in the numerical simulation. This confirms specific experimental observations according to which tumor cells are able to soften their nucleus and adjust their mechanical properties in order to facilitate extravasation [43].
In [12][13][14] a shape index (SI) parameter is used to quantify the nucleus self-deformation. SI is defined as follows: where S and l are the nucleus surface and perimeter, respectively. SI is equal to 1 for a perfect circle (i.e., no deformation), whereas it is equal to 0 for a straight line. In [12][13][14], SI is measured in the plan perpendicular to the micropillars. This is not possible in our model since we are in 2D. Nonetheless, we believe that it is even more interesting to measure SI in the sagittal plan to quantify the nucleus self-deformation, but also its penetration. We found that SI is equal to 0.32, 0.26 and 0.21 for S 1 , S 2 and S 3 , respectively indicating a significant nucleus deformation in particular when the nuclear lamina is ablated (i.e., S 3 ).

Conclusions
In this paper we have presented a 2D computational model to investigate the interactions between the cell and the micropillared substrate. The cell is initially suspended and gently meets the pillar due to the gravity force. Then, the contact force between the cell and the pillar is applied over a very thin layer.
The model is equipped with specific mechanical tools, namely the constitutive laws describing the cell components, their mechanical properties, and the fluid and solid phase mixture for each component. Our main objective has been establishing a correlation between such elements and the nucleus stress-strain state. To do so, a thorough sensibility study has been performed. More specifically, three series of simulations have been run: (i) the first one involves different materials models to describe the cell behaviour, i.e., Saint-Venant, Neo-Hookean, Mooney-Rivlin and Yeoh materials, (ii) the second one focuses on the balance between the fluid and solid phases of the lamina and nucleoplasm, (iii) the third one takes into account the variation of the mechanical parameters (i.e., the Young modulus) of the lamina and the nucleoplasm.
Through the large spectrum of combinations that we have provided, biologists may identify the one corresponding to specific cellular phenotypes or behaviors. Additionally, our model could inspire further experimental investigations to explore the interactions between the cell and the micropillared substrate.
We have been able to quantify the stress and the strain in the nucleus n . For the former, the Mises Cauchy stress r Mises has been assessed. For the latter, we have computed J 1 , J 2 and J 3 which provide the fibers elongation, the shear and length variations and the surface variations, respectively.
We have shown that the variation of the balance between the fluid and solid phase of the nucleus induces the maximum strains of the nucleus. In fact, we have found that varying c f ;l and c f ;np provides the highest values of J 1 ; J 2 and J 3 . More specifically, J 1 and J 2 are equal to 4.08 and 15.20 for c f ;l ¼ 0:7 and c f ;np ¼ 0:9 and such values are much higher than those found for different constitutive laws (2.68 and 5.244 for the Mooney Rivlin material) and for different values of E np and E l (2.44 and 4.27 for E np = 10 Pa and E l = 1000 Pa). For J 3 , the maximal values are closer although the highest value (1.31) is still found for the sensibility study on the fluid and solid phases balancing. This outcome confirms the importance of the phases balance and the need to correctly take them into account in the mechanical description of the nucleus. Furthermore, the amount of strain found could have an impact on the nucleus response and more particularly on the expression of strecth-sensitive proteins which could be responsible of cell motility in confined environments [40][41][42].
In terms of stress, even though the maximal values of the Mises Cauchy stress r Mises are very close for the three series of simulations, the highest one (1.47 Pa) is found when varying the Young moduli of the lamina and the nucleoplasm. More specifically, such a value is obtained when E np = 10 Pa and E l = 2000 Pa. This configuration could potentially induce the lamina rupture due to its high stiffness compared to the nucleoplasm. Such a phenomenon can change the lamina structure and therefore have critical consequences in disease like cancer or in genome stability [44][45][46].
Through some additional simulations, we have been able to compare our numerical results to the experimental observations on different cellular phenotypes. More specifically, we have shown that by increasing the simulation time, the nucleus undergoes significant self deformation and it acquires a 'peanut' shape by embracing the micropillar as it has been reported in [11][12][13][14]. Such a phenomenon is exacerbated when the nuclear lamina is ablated (i.e., reducing the Young modulus) and it confirms data according to which tumor cells are able to adjust their mechanical properties in order to invade healthy tissues [43].
To conclude, our model, which has been built only using mechanical tools, has drawn a large spectrum of scenarios to analyze and quantify the nucleus stress-strain state. It can help identifying the mechanical features responsible for specific nucleus responses and their impact on the global cell behavior. Although the interesting results, the present model could be improved in the different ways. First, a three-dimensional description of the system would allow to better evaluate the interactions between the cell and its surroundings and the nucleus role. Secondly, it could be interesting to quantify the pressure gradient in the cell in order to assess the proteins traffic flow to and from the nucleus. Finally, a precise description of the structure of the nuclear lamina could provide new information regarding its rupture and remodelling.
Funding Statement: The authors received no specific funding for this study.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.