A nonlinear elastic description of cell preferential orientations over a stretched substrate

The active response of cells to mechanical cues due to their interaction with the environment has been of increasing interest, since it is involved in many physiological phenomena, pathologies, and in tissue engineering. In particular, several experiments have shown that, if a substrate with overlying cells is cyclically stretched, they will reorient to reach a well-defined angle between their major axis and the main stretching direction. Recent experimental findings, also supported by a linear elastic model, indicated that the minimization of an elastic energy might drive this reorientation process. Motivated by the fact that a similar behaviour is observed even for high strains, in this paper we address the problem in the framework of finite elasticity, in order to study the presence of nonlinear effects. We find that, for a very large class of constitutive orthotropic models and with very general assumptions, there is a single linear relationship between a parameter describing the biaxial deformation and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos ^2\theta _{\mathrm{eq}}$$\end{document}cos2θeq, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{\mathrm{eq}}$$\end{document}θeq is the orientation angle of the cell, with the slope of the line depending on a specific combination of four parameters that characterize the nonlinear constitutive equation. We also study the effect of introducing a further dependence of the energy on the anisotropic invariants related to the square of the Cauchy–Green strain tensor. This leads to departures from the linear relationship mentioned above, that are again critically compared with experimental data.


Introduction
From the biological point of view, it is nowadays recognized that, in addition to chemical cues, cells actively respond to mechanical stimuli exerted on them by the surrounding environment. Even though the precise subcellular mechanisms governing this interaction still need to be understood, there is consistent evidence that the cytoskeleton has a central importance in the mechanical response (Hayakawa et al. 2001;Neidlinger-Wilke et al. 2001;Wang 2000;Wang et al. 2001). In particular, contractile actin stress fibers display the ability to develop reactive mechanical forces and to reorganize their structure in response to external changes: this mechanosensing is mediated by focal adhesions , that is, protein assemblies through which the actin internal structure is linked to the extracellular matrix (ECM). Focal adhesions provide both an anchorage to the substrate and a signal transmission, coupling the cell with the outside environment by detection of mechanical changes: they convert these stimuli into biochemical pathways, inducing a remodelling of the cytoskeletal structure.
These mechanical interactions between the cell and its neighbourhood are shown to play a fundamental role in several physiological situations, like cell motion, cell differentiation (Vlaikou et al. 2017) and tumour and tissue development. For instance, it is well known that cancer development is related to the loss of contact inhibition from its beginning (see Chaplain et al. 2006 and references therein) and that ECM stiffness and cell tensile stress influence both their proliferation and death (Butcher et al. 2009;Kass et al. 2007; Kumar and Weaver 2009). Incorrect response to mechanical cues is also related to many other diseases such as atrial fibrillation, intimal hyperplasia, scleroderma, diabetic nephropathy, glomerulosclerosis, enphysema, pulmonary, and liver fibrosis (Ingber 2009).
Stem cell fate and stem cell culturing are influenced by the mechanical interactions with the environment as well (Butcher et al. 2009;Discher et al. 2009;Guilak et al. 2009;Sun et al. 2012;Tulloch et al. 2011;Yan et al. 2013). This is the reason why cell cultures in vitro and tissue-engineered samples tend to reproduce the correct environmental conditions that real cells would live in (Costa et al. 2012). Consequently, just to mention two relevant examples, muscle cells and cardiac cells are often grown on substrata that are cyclically stretched, in order to mimic respectively muscle contraction and heartbeat (Kim et al. 1999;Laflamme and Murry 2011;Yoon et al. 2017). Now, when subject to a periodic deformation, several types of cells (not only those that have a bipolar morphology, such as fibroblasts, myofibroblasts, myocytes, airway smooth muscle cells, but also endothelial cells (Moretti et al. 2004)) show a peculiar response that proves their subtle sensitivity to mechanical prompts: when laying on a substratum they tend to re-orient themselves until they reach a stable configuration, with a well-defined angle between their polarization axis (along which the stress fibers become mainly aligned) and the direction of stretching. This reorientation process involves the cytoskeleton, with disruption and rebuilding of actin stress fibers, and the cell body, that follows the cytoskeletal reorganization with some delay (Hayakawa et al. 2001).
In the last two decades, cell reorientation following mechanical deformations has been investigated in an attempt to gain insight into this complicated phenomenon (Faust et al. 2011;Hayakawa et al. 2001;Jungbauer et al. 2008;Moretti et al. 2004;Morioka et al. 2011;Neidlinger-Wilke et al. 2001Wang et al. 2001Wang et al. , 2018. It is currently accepted that such a change in orientation is an active mechanism carried out by the cell (Wang et al. , 2018 and that mechanical strain is the main driving force of this process. However, experimental settings show a broad variety of behaviours: the vast majority of them proved that cells head their major axis towards a direction oblique or nearly perpendicular to the direction of greater stretch Hayakawa et al. 2001;Moretti et al. 2004;Neidlinger-Wilke et al. 2001;Wang et al. 2001), even though some authors reported different results (Bischofs and Schwarz 2003;Collinsworth et al. 2000). The final orientation angle reached by the cell is determined by several factors like the amplitude, the frequency Jungbauer et al. 2008), the biaxiality ratio of the applied cyclic deformation (Livne et al. 2014), and the substrate stiffness is sometimes thought to have a role as well (Obbink-Huizer et al. 2014;Tondon and Kaunas 2014;Tondon et al. 2012). For sake of completeness, we also mention that in 3D cell behaviour is influenced by the fact that they are embedded in a network of ECM fibers all around. So, under stretch, adhering to the fibers they tend to align more with the stretching direction (Asano et al. 2018;Chen et al. 2018;Eastwood et al. 1998).
From the modelling point of view, some attempts to explain this phenomenon have been done using linear elasticity to describe cell and substrate behaviour. Many works focused on looking for the directions of minimal strains or minimal stress (De et al. 2007(De et al. , 2008Faust et al. 2011;Livne et al. 2014;Wang 2000;Wang et al. 1995) as a preferred orientation for cellular placement under uniaxial cyclic strain conditions. However, in a recent work Livne et al. (2014) studied the response of cells on a substrate subject to biaxial extensions, finding a linear relation between cos 2 eq , where eq is the angle formed by the principal strain direction and the most elongated axis of the cell, and xx ∕( xx − yy ) , where xx and yy are the two infinitesimal principal strains of the biaxial test. The theoretical result correlated with the experimental data (Faust et al. 2011;Livne et al. 2014), suggesting that cells tend to minimize an elastic energy.
Nevertheless, it is worth observing that, in several experimental assays, mechanical tests were performed applying deformations for which using linear elasticity might be arguable, at least theoretically. As a matter of fact, for instance, in Livne et al. (2014) itself the maximum amplitude of the cyclic strain reached 24%, while in Faust et al. (2011) deformations up to 32% were applied to the specimen. In the latter work, the authors focused on the dependence of the orientation angle from the amplitude of the deformation, showing the existence of a stretch threshold necessary to induce a significant reorientation. Taken together, these results call for a study of the presence of nonlinear effects at high strains, in order to quantify their relevance in cellular mechanosensing and orientation dynamics. To our knowledge, only Lazopoulos and Pirentis (2007), Lazopoulos and Stamenović (2006) employed a finite elasticity framework to describe stress fibers reorganization in strained cells, although they considered only uniaxial substrate stretching and addressed the problem using a non-convex energy, giving an explanation based on the co-existence of phases.
The aim of this work is then to study the problem of cell reorientation in a nonlinear elasticity framework. The main goal is to understand why and to what extent the experimental results follow the same rule justified on the basis of linear elasticity also for large strains and are independent from the mechanical characteristics of the substratum and of the cells. We find that, considering the substratum with seeded cells as a finite-elastic orthotropic material, a very large class of elastic energies is minimized 1 3 by a relationship that can be considered as the nonlinear generalization of the one found in Livne et al. (2014). Specifically, introducing the parameter Λ ∶= x −1 x − y (where x and y are the principal stretches and x > y without loss of generality) that quantifies the biaxiality of the finite deformation, we put in evidence a linear relation between Λ and cos 2 eq , where eq is the equilibrium angle of the cell with respect to the x-axis. The slope of such a straight line turns out to be a combination of the coefficients that multiply the anisotropic invariants based on the right Cauchy-Green strain tensor ℂ and the polarization axis, which characterize the orthotropic constitutive model with quadratic dependency. The same relationship holds true for Fung-type materials as well, making the results valid for two energy forms often used in biomechanics.
Instead, a nonlinear dependence of the equilibrium from Λ comes out visibly if the elastic energy also depends on the anisotropic invariants related to the square of the Cauchy-Green strain tensor, i.e. ℂ 2 . In particular, we find that the presence of a non-vanishing coefficient in front of the term describing the response to stretch along the polarization axis gives rise to a departure from the linear behaviour that, for small values of the coefficient, is still compatible with experiments. Conversely, a non-vanishing coefficient in front of the term in charge of describing the response to stretch perpendicular to the polarization axis gives rise to results that look incompatible with experiments, suggesting this term is not present in the constitutive model. The dependence of the bifurcation point between the orientation perpendicular to the stretching direction and the oblique one is also examined for all cases, observing its disappearance above certain values of the energy coefficients.
In detail, the paper is organized as follows. In Sect. 2, we introduce the problem set-up and introduce the needed mathematical and mechanical background. In Sect. 3, we use the approach proposed in Saccomandi and Vianello (1997), Vianello (1996) to describe the stationary points for a general orthotropic constitutive equation. In Sect. 4, we focus on two very common types of energies employed in biological applications, namely a quadratic and Fung-type strain energy, showing that they behave identically as far as the stability problem is concerned. Afterwards, we calculate the specific equilibria and discuss their stability as a function of the involved material parameters. Section 5 is dedicated to the results in two significant cases, namely, in Sect. 5.1 we consider an energy which is independent of ℂ 2 , showing the existence of the linear relationship mentioned above, while in Sect. 5.2 we consider the effects of including the invariants that involve the square of the Cauchy-Green strain tensor. The final Sect. 6 summarizes the results and discusses some possible directions for future research.

Mechanical background and problem set-up
We consider a substratum provided with an ensemble of cells adherent to its surface in a subconfluent configuration, as in most experiments, which rules out the influence of cell-cell interaction. The system is then subject to a biaxial deformation, due to pulling and compression performed by an external device. As a consequence of this mechanical stimulus, following the deformation of the substrate cells orient on average along a certain direction, which can be identified through a unit vector in the first quadrant of a reference system with the axes along the directions of the biaxial deformation (see Fig. 1). If we assume that the system composed by the substrate and the overlying cells behaves as an elastic continuum, the deformation will induce a storage of energy into the body that depends on the orientation . We want here to study how the elastic strain energy of the system subject to a biaxial stretch depends on the average direction of the stress fibers .
For this purpose, we consider a general elastic energy density U for an orthotropic material depending on the deformation gradient through the invariants of the right Cauchy-Green deformation tensor ℂ = T where ⟂ is the direction perpendicular to in the plane of the substratum, which is also the plane containing the principal stretch. Notice that in (1) one could also expect a dependence on ⟂ ⋅ ℂ 2 , but it is omitted since it depends on the other invariants (see, for instance, Ogden 2003).
This kind of energy is commonly used to describe the mechanical behaviour of anisotropic materials that show two preferential directions (Holzapfel and Gasser 2001;Ogden 2003). For instance, it is employed for fiber-reinforced materials, in which there are, say, two orthogonal bundles of fibers that influence the mechanical response of the body (Melnik and Goriely 2013), or for blood vessel walls that present fiber bundles in different directions. In our case, the system displays a natural anisotropy due to the presence in the cells of aligned stress fibers (SF) and actin filament structures that are cross-linked by several types of proteins, like fascin, fimbrin, -actinin, filamin, ARP2-3 (Civelekoglu-Scholey et al. 2005;Lu et al. 2008;Xu et al. 2016) (see Fig. 1). Since SF are two to even ten times stiffer than the lateral actin network (1) U = U(I 1 , I 2 , I 3 , I 4 , I 5 , I 6 , I 7 , I 8 ), (2) (Lu et al. 2008;Mathur et al. 2007), they induce bidirectional anisotropy in the mechanical response, justifying the general assumption (1) to consider the system as orthotropic.
We also observe that a cell does not have a real polarization given by a head and a tail (Wang et al. 1995). So, configurations with cells aligned along and − are geometrically indistinguishable and therefore energetically equivalent. A similar thing holds true if we replace ⟂ with − ⟂ . This implies that the elastic energy should be the same under the related symmetry transformations. All the invariants mentioned in (2) achieve the same values under the symmetry transformations above, except for I 8 that changes sign, i.e. given ℂ , I 8 (− , ⟂ , ℂ) = −I 8 ( , ⟂ , ℂ) and I 8 ( , − ⟂ , ℂ) = −I 8 ( , ⟂ , ℂ) . Therefore, to satisfy the symmetry requirements, U must be an even function of I 8 .
If we orient the x-and y-axes of the reference frame along the principal stretching directions, the right Cauchy-Green tensor is diagonal and its eigenvalues are the principal stretches, i.e. Fig. 1 Sketch of experimental set-up and consequent inner structure of a typical cell with x > y and x > 1 . Actually, in the case of equi-biaxial deformation x = y , cells do not show a preferential orientation . Once the deformation is fixed, our goal is to study which orientations for the stress fibers correspond to minima of the elastic energy.

Stationary points by the coaxiality approach
The general problem of finding the critical points of an hyperelastic anisotropic strain energy for a body subject to a rotation and a deformation has been studied in Saccomandi and Vianello (1997), Vianello (1996). In particular, it has been shown that the critical points of the energy are achieved for those rotations that make the stress and strain tensors coaxial. Since two symmetric tensors are coaxial if and only if they commute (Vianello 1996), we have to find the rotations ℚ about the z-axis such that where ℂ * = ℚℂℚ T and * = (ℂ * ) is the second Piola-Kirchhoff stress tensor corresponding to the deformation ℂ * .
To explicitly write * in our case, we define the structural tensors and recall that Therefore, the second Piola-Kirchhoff stress tensor reads where we have defined Then, using (6) If we focus on the last term on the left-hand side, elementary tensor algebra allows to rewrite it as with while the operator Sym takes the symmetric part of its tensorial argument. A direct substitution into (9) leads to or equivalently where It is straightforward to observe that the conclusion is not influenced by the classical first three invariants I 1 , I 2 and I 3 , which represent the isotropic response of the material and do not change with rotations. Hence, since we are interested in the anisotropic behaviour of the system as a consequence of cell orientation, their dependence will be dropped in the following discussion, though one should remember that these (8) terms might appear in any coefficient or as an additional term in the energy.
Equation (11) states that stationary points are identified by the rotations ℚ about the z-axis which make the tensors ℂ and ̂ defined in (12) commute. As pointed out in Vianello (1996), there are always at least two solutions, which can be easily identified in this case. In fact, if ℚ is such that ℚ is along one of the principal directions, then, on the one hand, * = ℚ T ℚ is diagonal and therefore commutes with ℂ , which is diagonal by definition. On the other hand, we have that I 8 ( , ⟂ , ℂ * ) = I 8 (ℚ T , ℚ T ⟂ , ℂ) which vanishes. Consequently, s 8 (I 8 ) is null as well because it is an odd function of I 8 , because of the symmetry requirements on U which is assumed to be even in I 8 .
We conclude that configurations with cells oriented along the principal stretching direction or perpendicularly to it always correspond to stationary points of the elastic energy.
However, one might have further equilibria for other possible rotations satisfying (11), which will depend in general on the specific energy functional considered. In the following we will show that, for a large class of elastic energies, there might be two symmetric equilibria and we will study the stability of all configurations.

Stability conditions for a quadratic-like energy
If we use the reference frame with the axes along the principal stretching directions, so that ℂ is diagonal, it is convenient to identify the cell major axis through the angle it forms with the x-axis, so that = (cos , sin , 0) and ⟂ = (− sin , cos , 0) . This angle is univocally associated to a rotation ℚ in the framework of the previous Section and in the following will be used as our main variable. With this choice, the invariants in (2) read and can be compactly rewritten as where (13) I 4 = x cos 2 + y sin 2 = ( x − y ) cos 2 + y ,

The elastic energy
In order to investigate the existence of other stationary configurations in addition to the trivial ones, and to study their stability, we need now to specialize the elastic energy, trying to keep it as general as possible. We consider then the class of elastic energies that can be written as a homogeneous second-order polynomial in the variables Î i ∶= I i − 1 , for i = 4, 5, 6, 7 , and I 8 , plus a term related to the isotropic response: where ∶= Î 4 ,Î 5 ,Î 6 ,Î 7 , I 8 and is the symmetric matrix of coefficients (that might possibly depend on the isotropic invariants), while V is the purely isotropic contribution that depends on (I 1 , I 2 , I 3 ).
We observe that the following analysis can be straightforwardly repeated for a Fung-type energy that is often used in biomechanical applications (Fung 1967), giving rise to the same results. In fact, the stationary points of U F coincide with the ones of U . Moreover, their stability character may be identified by the second derivative of U and then is the same as well. Therefore, the results obtained for a quadratic energy also hold for an exponentiallike energy, amplifying the validity of our conclusions.
We also remark that, in order to slightly reduce the terms influencing the stability analysis, we do not consider here possible isotropic-anisotropic couplings in the energy (15), that is, we exclude terms like I hÎl , where h ∈ {1, 2, 3} and l ∈ {4, 5, 6, 7, 8} . This can be easily done and their role will be analysed in a future work.
For future convenience, it is useful to denote by k ij for i, j = 4, … , 8 the coefficients of the matrix (e.g., k 44 is the coefficient in the top left corner of the matrix). In particular, the coefficient k 44 is related to the stiffness along the polarization axis of the cell and k 66 to the one along the direction orthogonal to the cell major axis. Considering that stress fibers are mainly aligned to the polarization direction, coherently with Butler et al. (2002) we will assume that We also point out that the coefficient k 88 is related to the response to shear, i.e. at the microscopic level to the resistance of changing angle among actin fibers, also involving the action of actin cross-linking proteins, such as filamin, Rho/ Rac GTPases (Civelekoglu-Scholey et al. 2005;Wang et al. 2001) and Arp2/3 as well (Goley and Welch 2006;Rouiller et al. 2008) (see Fig. 1).
We finally observe that, because of the symmetry conditions related to a switching of the orientation of the axis, the coefficients k i8 (and k 8i ), i = 4, 5, 6, 7 must vanish. In fact, they multiply the cross terms Î i I 8 that change sign under these symmetry transformations. Then, their inclusion would break the symmetry of the energy, leading to a biologically unfeasible situation.

Stationary configurations
Taking into account the consequences of symmetry on the coefficients of matrix in (15), the stationary configurations are then identified by the solutions of that can be explicitly written as Recalling the form of the coefficients (14), the previous equation reads As already discussed at the end of the previous Section, the configurations with = 0 and = ∕2 (which in the following will sometimes be referred to as parallel and perpendicular orientation, respectively) are always stationary. However, there might be additional equilibria satisfying or where that can also be explicited as k ij a i a j , and The nontrivial solution eq , also called oblique orientation, exists if and only if or equivalently if

Stability
As far as stability is concerned, we need to evaluate the second derivative of the energy Imposing the positivity of U ′′ then leads to the following stability conditions: where in the last condition we have used the fact that U �� ( eq ) > 0 requires Σ 1 − k 88 > 0 to simplify the existence condition (19). In the following, we will sometimes refer to the preferential orientations of Eq. (20) as equilibrium angles. However, we remark that they are not in general solutions of the elastic problem, but only stationary points of the elastic energy.

Bifurcation results
Starting from the computations performed above, in this Section by means of a bifurcation analysis we discuss the preferential orientations of cells on a stretched substrate, for k ij a j (a i cos 2 + b i )(cos 2 − sin 2 ) + k 88 a 2 4 (cos 2 − sin 2 ) 2 − 4 cos 2 sin 2 . (20) an elastic energy given by (15) or (16). We firstly discuss the case in which the strain energy does not depend on the invariants related to ℂ 2 , and then we analyse the effects of their correction, since they introduce a significant change in the predicted behaviour. In particular, we consider three types of finite deformation: in the first one we fix the stretch in the y-direction y to a certain value and z = 1 , while letting x vary. In the second one we keep x y = 1 , corresponding to an isochoric planar deformation if z = 1 . Finally, we set y = 1∕ √ x , equivalent to a pure isochoric deformation that involves also the z-direction. Although we do not have any experimental information on z , it does not affect the discussion on equilibria or on their stability character.

Independence from ℂ 2
Let us firstly consider the case in which the elastic energy is independent from terms containing ℂ 2 , that is, from the anisotropic invariants I 5 and I 7 . This case can be easily obtained from the previous computations setting k 45 = k 47 = k 55 = k 56 = k 57 = k 67 = k 77 = 0 . Doing so, we can get simplified forms for the terms and Hence, in addition to the trivial equilibria, Eq. (17) simplifies to putting in clear evidence a linear relationship between cos 2 eq and the parameter which compares the elongation along x with respect to the sum of elongation and contraction along y.
The slope of the straight line (21)  It is also worth to observe that the line always passes through the point corresponding to cos 2 eq = 1∕2 (i.e. eq = ∕4 ) when Λ = 1∕2 , for any values of the elastic parameters. This is explained by the fact that, for a deformation satisfying Λ = 1∕2 (or equivalently x − 1 = 1 − y ), the minimum of the energy coincides with the direction of minimal strain. Indeed, the latter is given by the angle such that I 4 = 1 , that is, The black lines refer to the case | | > 1 and the red lines to | | < 1 , while dashed lines represent unstable configurations and full lines stable ones. In (b) full lines indicating the stability of = 0 till the bifurcation point and of = 2 for < −1 are not drawn. The area with Λ < 0 is not shown because of the pulling characteristics of the experiments Therefore, when Λ = 1∕2 , we have ̃ = ∕4 which coincides with eq obtained using Eq. (21). This is not true in general, since ̃ satisfies differently from what is found through energy minimization in Eq. (21) except for the case = 1 , which does not fit the experimental data (Faust et al. 2011;Livne et al. 2014). In fact, as pointed out in Livne et al. (2014), choosing the minimal strain direction as the preferential one for cell orientation does not allow to describe the experimental observations, while an energetic approach does. The only case in which the minimal strain direction and the minimal energy direction coincide is precisely when Λ = 1∕2 . These observations also allow us to justify even more the need of choosing an orthotropic constitutive model instead of a purely transversely isotropic one. In fact, the elastic strain energy for a transversely isotropic medium only depends on I 4 and I 5 in addition to the three isotropic invariants, and then only on I 4 when the dependence on ℂ 2 is not considered as done in this Section. However, in such a case one obtains = 1 and so minimizing a transversely isotropic energy would be equivalent to choosing the minimal strain direction again, which as just discussed does not seem accurate.
We also remark that Λ = 1 corresponds to clamping the specimen, so that y = 1 , while Λ > 1 corresponds to stretching also along y, still keeping y < x . On the other hand, values of Λ < 1 2 correspond to x − 1 < 1 − y , i.e. compressions along y are stronger than elongations along x, which is not done in experiments reported in the literature.
Finally, we observe the following cases related to isochoric deformations, which will be examined in detail later: with Λ → 1 for very large x and the lower extremum of the interval achieved in the limit of no stretching. Of course, for these types of finite deformation, the case Λ > 1 is unfeasible.
In terms of and Λ the existence condition for the nontrivial equilibrium (18) writes as

Case ˛> 0
When > 0 , the stability conditions (20) become that is, the non-trivial equilibrium position is stable whenever it exists.
Recalling that k 44 > k 66 and focusing for the moment only on the case > 0 , there are two relevant subcases to discuss, depending on the following relationships between parameters:

Referring to Fig. 2a, in case (ii)-represented by red lineswe have two supercritical bifurcation points at
The configuration with orientation perpendicular to the stretching direction (i.e. eq = ∕2 , cos 2 eq = 0 ) is stable if Λ > 1+ 2 . When Λ is decreased below this value, the polarization axis of the cell tends instead to orient obliquely, finally becoming completely aligned to the stretching direction if Λ < 1− 2 . However, in case (i), represented by the black line in Fig. 2a, this value is negative and so it cannot be physically achieved in the usual experimental set-up. In this situation, cells will never orient themselves along the stretching direction.
In order to compare this behaviour with previous linear elasticity results, we observe that if we define and r such that x = 1 + and y = 1 − r , we have that Λ = 1 1+r . So, one recovers the linear relationship between cos 2 and 1 1+r discussed in Livne et al. (2014). Actually, since = 0.794 ± 0.08 < 1 , the situation observed in their experiments seems to correspond to case (ii). In this case, the two bifurcation points fall in the interval [0, 1]. We recall that Λ = 1 corresponds to y = 1 , that is no stretching in the y-direction, i.e. a clamped condition.
The analysis done here shows, in particular, that any model of type (15) or of Fung type (16) independent of I 5 and I 7 (i.e. the invariants depending on ℂ 2 ) with the following relation among the coefficients is able to fit the parameters, in a way that is independent of the magnitude of the applied strain, even outside the range of validity of linear elasticity. This explains why the experimental behaviour shown in Faust et al. (2011), Livne et al. (2014) seems to be independent or nearly independent of the magnitude of the applied strain.

Case ˛< 0
For the sake of completeness, we analyse the case in which < 0 , which can occur, for instance, if k 88 is much larger than the other parameters. In this case, noticing that 1+ 2 < 1− 2 , the stability conditions are the following: More precisely, also in this case there are two distinct situations to be considered: Referring to Fig. 2b, in case (iii)-corresponding to red lines-one has two admissible subcritical bifurcation points with coexistence of two critical equilibria if Λ ∈ 1+ 2 , 1− 2 , corresponding to the parallel and perpendicular orientation, and only one of them stable outside this range. In case iv) the equilibrium = ∕2 is always stable while = 0 is stable only if Λ < 1− 2 . In all cases, if < 0 the oblique orientation is unstable.

ℂ 2 −correction
We now analyse the correction that will be introduced if a dependence of the energy on the invariants depending on ℂ 2 (i.e. I 5 and I 7 ) is allowed. For this purpose, it is convenient to rewrite with (23) k 44 + k 66 − 2k 46 − k 88 k 44 − k 66 ≈ 0.794 and where After computation, we can explicit the nontrivial equilibrium (17) as We notice that the terms A 1 , A 2 and B 2 are the corrections related to the presence of ℂ 2 -dependent terms through the invariants I 5 and I 7 . Indeed, when the elastic energy does not include their contribution, we recover (21), i.e. the linear relationship between cos 2 and Λ . However, since the correction coefficients depend on x and y , Eq. (25) does not represent a straight line anymore. It is worth remarking that the last term in (25) represents a shift from the equilibrium angle = ∕4 that looks close to the one observed when Λ = 1∕2.

Small deformation limit
We observe that, in the limit of small deformation, the zeroth order approximations of the corrections become and in Eq. (25) the last term vanishes, so that the approximating curve is again a straight line: (25) such that eq = ∕4 when Λ = 1∕2 (that corresponds to x − 1 ≈ 1 − y ≪ 1 ). Hence, in the linear limit, the only difference from the case analysed in the previous Section is that more coefficients contribute to the identification of the slope, or viceversa, given some experimental data it is hard to distinguish which coefficients contribute to the slope of the line. This is not surprising: in the linear limit, for instance, the contribution to the energy of I 4 is indistinguishable from that of I 5 , and the dependence on I 6 merges with the one on I 7 .
In fact, Eq. (26) can be rewritten in the same form as (21) with the following formal substitutions: In the following, in order to evaluate the influence of ℂ 2 corrections, we change one parameter at a time while keeping the slope of the straight line in (26) constant, in order to start from the same linear dependence. For instance, when k 55 > 0 , then, recalling (27), k 44 is decreased accordingly, so that the value of K ‖ is mantained constant. In particular, with this idea in mind, we fix the coefficients K ‖ = 0.4 , K ‖⟂ = 0.1 , K ⟂ = 0.1 and k 88 = 0.0618 in order to match the experimental fitting value reported in Eq. (23).

Effect of k 55
In Fig. 3 we focus on the effect of a non-vanishing k 55 , keeping to zero the other coefficients involving the indices 5 and 7. As already stated above, to keep the same value of K ‖ and therefore the same linear limit given by (26), as we vary k 55 we accordingly change k 44 = K ‖ − 4k 55 . To observe how the non-trivial equilibrium is changed, we focus on the three types of deformation defined at the beginning of Sect. 5. One can observe that in all cases changing k 55 leads to a departure of the equilibrium curve from a straight line, becoming convex. In addition, keeping y fixed and changing Λ as in Fig. 3a, one can appreciate a decrease in the value of cos 2 eq , which means an increase in the equilibrium angle, for Λ = 1∕2. x . For all the three types of deformation, we see that the presence of k 55 leads to a departure from the linear relation between cos 2 eq and Λ . In d the bifurcation point Λ b for which eq = ∕2 is shown as a function of k 55 : the plot highlights that, for sufficiently high values of k 55 (for instance k 55 > 0.085 for y = 0.8 ), the bifurcation point disappears, since values of Λ b > 1 are not admissible for the deformations we consider On the other hand, when x y = 1 as in Fig. 3b, eq stays fixed at ∕4 if Λ = 1∕2 . In fact, in the latter case A 1 = A 2 = B 2 = 4k 55 and the numerator of the last term in (25) vanishes for Λ = 1∕2 . This does not occur in the former case. At the other extreme, when Λ → 1 , then the denominator of the second term in (25) grows indefinitely, while the last term tends to 1/2. So, the curve given by (25) tends back to 0 when Λ → 1 (not visible in the figure because the curve becomes negative).
Finally, in the third case y = √ x shown in Fig. 3c we have a minimal possible value of Λ = 2∕3 , corresponding to an equilibrium angle such that cos 2 eq = 0.29 , i.e. eq ≈ 57 • , that is sometimes reported in some discussions of experiments (see, for instance, Wang 2000; Wang et al. 1995).
As regards the bifurcation point between the equilibrium branch = ∕2 and the one given by (25), it is implicitly defined by since the coefficients A 2 and B 2 also depend on Λ b . Therefore, while we have seen that the introduction of nonlinear terms depending on I 4 , I 6 and I 8 does not modify the bifurcation obtained in the linearized theory, the dependence of the energy on the invariants related to ℂ 2 entails an identifiable shift of the bifurcation point between the perpendicular orientation and the oblique one that may also disappear for large values of k 55 , as shown in Fig. 3d.
If we assume k 55 ≠ 0 , the shift of the bifurcation point for eq = ∕2 can be evaluated through This means that, for given coefficients but k 55 , (31) implicitly defines Λ b in terms of k 55 that, actually, can be made explicit by writing We notice that, in the equation above, the term (1 + y )( x + y ) is a function of Λ b , given, for instance, by and by .
If y = 1∕ √ x we have that which corresponds to These relations can be exploited to plot the variation of the bifurcation point as a function of k 55 , as shown in Fig. 3d. The bifurcation point Λ b approaches 1 for increasing values of k 55 , eventually leading to a disappearance of the bifurcation for some critical values of the parameter. To be more specific, if y = 0.8 , for k 55 > 0.085 the branch relative to the oblique equilibrium and = ∕2 do not cross for physically admissible values of Λ and therefore the perpendicular orientation is always unstable. In particular, we observe that in the case of fixed y the threshold value is higher, while for the other two types of deformation it is the same and amounts to k 55 = k 88 ∕2 . This follows immediately from (32) and (33) for y = 1∕ x , while in the case y = 1∕ √ x it is enough to observe that Substituting into (32) and recalling that Λ b → 1 means x → +∞ immediately leads to k 55 = k 88 ∕2 , coherently with Fig. 3d.
In order to evaluate the relevance of the nonlinear correction, we performed a fitting of experimental data extracted from Livne et al. (2014) where, using a biaxial experiment, the stretches in the two directions are actually controlled. Although we do not know the exact value of y , we observe that if y = 1 , then Λ would be identically equal to 1. So, we fix it to be y = 0.8 . Then, focusing on Fig. 4, we explored the possibility of a better fitting using nonlinear elasticity and ℂ 2 correction. Actually, as already observed, data are already fitted quite well by the straight line. However, using a nonlinear regression estimation we find that a small value of k 55 = 0.008 gives an even better fitting for data in Livne et al. (2014). We also find that, if we increase the fixed value of y , a higher value of the coefficient k 55 is needed to fit the data. For instance, when we take y = 0.98 we find a best fitting value of k 55 = 0.08 , an order of magnitude greater than the one obtained for y = 0.8 . This is due to the fact that, when y ≈ 1 , a small stretch x in the x-direction is sufficient to span all the ( admissible values of Λ . Consequently, nonlinear effects become less relevant and to fit the nonlinear model we need to take very high values of the related coefficients. In Fig. 5, we present the fitting using also the data by Faust et al. (2011): we explored the fixed y case since, even if the Poisson ratio of the specimen is taken to be = 1∕2 , the actual reported biaxiality ratio is different. In particular, considering the experimental settings with 32% and 31.7% strain, one finds a value of y = 0.952 and y = 0.91 , respectively, since the biaxiality ratio is Λ = 0.87 and Λ = 0.78 in the two cases. In Fig. 5a, we then tried to fit the data (Faust et al. 2011;Livne et al. 2014) simultaneously to assess if a nonlinear correction could better explain the experimental observations in different settings. It is found that, for a fixed value of y = 0.952 , taking k 55 = 0.04 gives a slightly better fitting than the straight line approximation, while nonlinear regression performed only on the data in Livne et al. (2014) returns a value of k 55 = 0.03 . To highlight this difference, in Fig. 5b we plot the deviation from the linear approximation, that is c Equilibrium orientation as a function of k 55 for x = 1.32 and y = 0.952 , related to experimental actin orientations obtained by Faust et al. (2011) for a 32% stretch (represented by the circle). As shown, a value of k 55 = 0.008 is able to capture the experimental orientation precisely 1 3 as a function of Λ . Since most of the experimental data fall below the straight line, i.e. below 0 in Fig. 5b, a convex curve obtained with the introduction of k 55 is able to better approximate the observed behaviour, even if the difference is of the order of 10 −2 . Finally, we chose one of the experimental actin angle values obtained in Faust et al. (2011) for a 32% stretch and fixed x and y in order to have the same value of Λ used in the experiment. Doing so, we tried to find a value of k 55 able to capture this experimental point: as shown in Fig. 5c, a value of k 55 = 0.088 precisely fits the orientation angle for such a fixed deformation.
We conclude our analysis by making a comparison with the transversely isotropic case, which we already discussed as inadequate to fit the experimental data. This is confirmed by the curves reported in Fig. 6. Indeed, if we consider a transversely isotropic energy that depends only on five invariants, recalling (25) we get: In Fig. 6a we plot the relationship (37) for different values of k 55 and k 45 = 0 , while in Fig. 6b the case k 45 ≠ 0 is shown. It is clearly observed that, in both cases, the transversely isotropic model provides a fitting which is not satisfactory compared to the orthotropic one reported in Fig. 4. To have a better insight, in Fig. 6c we show a direct comparison between the best fitting curves in the transversely isotropic and orthotropic case: in the latter, there is a significant improvement in the fitting of experimental data.

Effect of k 77
We turn now the attention to the effect of a non-vanishing k 77 , keeping K ⟂ fixed, and perform a similar reasoning as we did for k 55 . Again, as shown in Fig. 7b, if x y = 1 then = ∕4 when Λ = 1∕2 for any value of k 77 . In fact, in this case, A 1 = 4k 77 and B 2 = −4k 77 , while A 2 = 0 . So, the numerator in the last term of (25) vanishes. This does not occur if y is kept fixed, as shown in Fig. 7a. In fact, the value of eq slightly decreases for increasing values of k 77 . A more dramatic effect occurs if k 77 ≠ 0 when Λ → 1 , because as before the second term in (25) tends to zero, but on the contrary of the previous case the last term tends to −1∕2 . This implies that the curve in (25) tends to 1 when Λ → 1 and, as a consequence, there are two bifurcation for different values of k 55 is shown, while in b the effect of k 45 is investigated. In c a direct comparison between the best fitting curves for the transversely isotropic and orthotropic case is provided, showing that the latter is more accurate values for the equilibrium configuration eq = ∕2 that is stable for values of Λ in between them. However, the two bifurcation points soon disappear. Specifically, as shown in Fig. 7d, this happens for k 77 > 5.4 ⋅ 10 −4 , if y = 0.8 fixed.
For such values, the equilibrium eq = ∕2 is always unstable, while the branch represented in Fig. 7 is always stable. This does not seem to correspond to what is observed in experiments, suggesting that k 77 = 0 or that it must have a very small value in the constitutive model.

Effect of k 45
Finally, we focus on the effect of coupling terms other than k 46 , that is already present in the analysis of the previous Section, as it involves the right Cauchy-Green tensor ℂ . Having ruled out the importance of a non-vanishing k 77 we do not report the plots involving k 47 , k 57 and k 67 and focus instead on the remaining mixing parameter k 45 . The results of its introduction in the model are shown in Fig. 8. In particular, we fixed a value of k 55 = 0.01 and plotted the equilibrium angle as a function of Λ for the three different deformations mentioned above. We observe an effect on the equilibrium curve similar to the one encountered increasing k 55 . The departure from the straight line is slighter if y = 0.8 , and greater in the other two cases. Focusing on Fig. 8d, as in the case k 55 ≠ 0 , the bifurcation value above which the configuration eq = ∕2 is stable increases with k 45 and disappears for higher values of this parameter because the curve corresponding to the non-trivial configuration decreases to 1 (without a minimum). As a last observation, comparing the theoretical results with experimental data shown in Fig. 8a it does not seem that the introduction of the mixing term k 45 is relevant for a better fit, as confirmed by nonlinear regression.

Discussion
Stimulated by the need of understanding why the preferential configurations achieved by an ensemble of cells over a substratum under stretch were apparently independent of the mechanical characteristic of the substratum itself, of the stretch magnitude, and to a certain extent of the cell type, we investigated the phenomenon from a mechanical point of view. We then worked in the framework of nonlinear elasticity with a very general class of constitutive equations, often used in biomechanical applications, to describe the cell and the stretched substrate as an orthotropic material undergoing finite deformations. Fig. 7 Analysis of non-trivial equilibrium position for k 77 ≠ 0 and k 66 = 0.1 − 4k 77 to keep the same linear limit. The other non-vanishing parameters are k 44 = 0.4 , k 46 = 0.1 , and k 88 = 0.0618 . In a y = 0.8 , while in b y = 1∕ x and in c For all the three types of deformation, one can see that introducing k 77 provokes a significant difference from the linear behaviour not observed in experiments, suggesting that such a parameter is not present in the constitutive model. In d the bifurcation points Λ b for which eq = ∕2 are shown as a function of k 77 . The region enclosed by the curve identifies the values of Λ for which eq = ∕2 is stable. For sufficiently high values of k 77 the bifurcation points disappear and the configuration eq = ∕2 is always unstable Coherently with experiments and previous theoretical works that use linear elasticity, our model predicts that strained cells will reorient until their major axis forms a well-defined angle with the direction of stretching, depending on the biaxiality of the deformation quantified by the parameter Λ . We found that the bifurcation diagram only depends on a particular combination of the coefficients k 44 , k 66 , k 88 (in addition to the mixed one k 46 ) appearing in front of the terms with the anisotropic invariants I 4 , I 6 , and I 8 , that is, the ones determining the influence of ⋅ ℂ , ⟂ ⋅ ℂ ⟂ , and ⟂ ⋅ ℂ , where is the direction of cell polarization and ⟂ the perpendicular one.
Hence, the results obtained using linear elasticity extend their validity also in the nonlinear regime and for a wide class of orthotropic models, which explains why experimental results with a stretch amplitude outside the linear elastic regime can be still well described by the linearized theory. In particular, our predictions confirm that there exists a linear relationship between Λ and cos 2 eq , where eq is the angle formed by the aligned stress fibers of the cell and the main stretch direction, with slope that from experiments is given by Therefore, we showed that any quadratic-like or Fung-type constitutive model independent of the invariants related to ℂ 2 (i.e. I 5 and I 7 ) can fit the experimental data, provided that the coefficients satisfy the relation (38). The only difference from the theoretical point of view is that, using finite elasticity, one has more coefficients contributing to the slope of the bifurcation line. Besides, our analysis is in principle able to describe the final reorientation angle even in the cases < 0 and > 1 , which might be of interest in other experimental settings.
We also proved that for Λ > 1+ 2 ≈ 0.897 the only stable equilibrium is the one in which cells are aligned perpendicularly to the main stretching direction.
Then, we introduced the corrections due to the invariants related to ℂ 2 and analysed their influence. It is found that the presence of the invariant I 5 involving ⋅ ℂ 2 has an effect which is still compatible with experimental data, for small values of the pertaining coefficient k 55 and for all the three types of deformation considered. In this case, however, the dependence discussed above is no longer linear (38) k 44 − k 66 k 44 + k 66 − 2k 46 − k 88 ≈ 1.26 . Fig. 8 Analysis of nontrivial equilibrium position for k 45 ≠ 0 , considering k 55 = 0.1 and k 44 = 0.4 − 4k 55 − 4k 45 to keep the same linear limit. The other non-vanishing parameters are k 46 = 0.1 and k 88 = 0.0618 . In a we set y = 0.8 , while in b y = 1∕ x and in c y = 1∕ √ x . Moreover, in a we reported the experimental data obtained in Livne et al. (2014): the introduction of the mixing parameter k 45 does not improve significantly the fitting of such data. In (d) the bifurcation point Λ b for which eq = ∕2 is shown as a function of k 45 , when k 55 = 0.01 is kept fixed. For all the types of deformation considered, the bifurcation disappears when k 45 crosses a critical value, which is the same for y = 1∕ x and and the effect due to finite elasticity comes out clearly. This allows us to speculate that there might be a slight further nonlinear impact in the phenomenon, since the introduction of k 55 provided a slightly better fit to experimental data. The presence of such a parameter also influences the existence of a bifurcation between the perpendicular and the oblique orientation, which disappears for sufficiently high values of k 55 . However, it is worth to remark that the data which are present in the literature seem very different: there is then a degree of uncertainty which does not allow to precisely conclude on the relevance of such a dependence. More experiments would be needed to assess the exact influence of this nonlinear parameter, progressively varying the stretching amplitude. From our results, it seems that performing experiments at values of Λ close to one, i.e. nearly clamped experiments in the y-direction, would better put in evidence the effects of I 5 , as well as focusing on the moment when cells start departing from a perpendicular orientation with respect to the stretching direction, i.e. the bifurcation point occurring for high values of Λ. Instead, it is found that a dependence of the strain energy on the invariant I 7 involving ⟂ ⋅ ℂ 2 ⟂ yields results that look incompatible with experiments. In particular, we showed that the presence of the coefficient k 77 leads to a consistent difference from the linear behaviour especially for high values of Λ and to the disappearance of the bifurcation between the perpendicular and oblique orientation. Since these predictions are not confirmed by experimental data, we speculate that the constitutive model should not include its contribution. Finally, we showed that the influence of a mixing parameter related to ℂ 2 like k 45 , which couples the effect of the invariants I 4 and I 5 , is negligible and does not provide a better fitting to experimental data.
An aspect that is not covered by the present modelling is the dependence of the behaviour on the frequency of the applied cyclic stretch. In order to include that, one should take into account the viscoelastic properties of the material and in particular the characteristic response times of the cells with respect to the period of imposed oscillations. Some works in this direction have been done in Kong et al. (2008), Qian et al. (2013), who considered the viscoelasticity of actin stress fibers, mostly from a microscopic viewpoint, while in an ongoing work we are considering linear viscoelasticity regarding the system as a continuum.
More importantly, during the experimental tests there are remodelling effects related to the internal reorganization of the cytoskeleton and of the stress fibers. Introducing remodelling into the continuum structure of the model as done in Ciambella and Nardinocchi (2019), Ciambella and Nardinocchi (2020) could allow to explore this effect, maybe improving the description of such a complex behaviour of the cells.