A Statistical Mechanics Approach to Describe Cell Reorientation Under Stretch

Experiments show that when a monolayer of cells cultured on an elastic substratum is subject to a cyclic stretch, cells tend to reorient either perpendicularly or at an oblique angle with respect to the main stretching direction. Due to stochastic effects, however, the distribution of angles achieved by the cells is broader and, experimentally, histograms over the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[0^\circ , 90^\circ ]$$\end{document}[0∘,90∘] are usually reported. Here we will determine the evolution and the stationary state of probability density functions describing the statistical distribution of the orientations of the cells using Fokker–Planck equations derived from microscopic rules for describing the reorientation process of the cell. As a first attempt, we shall use a stochastic differential equation related to a very general elastic energy that the cell tries to minimize and, we will show that the results of the time integration and of the stationary state of the related forward Fokker–Planck equation compare very well with experimental results obtained by different researchers. Then, in order to model more accurately the microscopic process of cell reorientation and to shed light on the mechanisms performed by cells that are subject to cyclic stretch, we consider discrete in time random processes that allow to recover Fokker–Planck equations through classical tools of kinetic theory. In particular, we shall introduce a model of reorientation as a function of the rotation angle as a result of an optimal control problem. Also in this latter case the results match very well with experiments.


Introduction
In the 80's the study of cardiovascular diseases led to the need of understanding the behaviour of cells of the heart and of the arterial walls subject to periodic deformations due to pulsatile heart contractions and consequent blood flow (Buck 1979(Buck , 1980. In order to mimick this environment, many authors seeded cells on a substratum that was stretched periodically (see, for instance, the recent review (Giverso et al. 2022) and references therein). It was generally found that for sufficiently high stretching frequencies (see Greiner et al. 2015;Hsu et al. 2009;Jungbauer et al. 2008;Lee et al. 2010;Tondon and Kaunas 2014) and amplitudes (see Boccafoschi et al. 2007;Dartsch et al. 1986;Kaunas et al. 2005;Mao et al. 2021;Morita et al. 2013), cells internally develop stress fibers that link to the substratum via focal adhesions and confer anisotropic characteristics to the cell (see Fig. 1). Such stress fibers are, at the equilibrium, mainly aligned perpendicularly to the main stretching direction or at oblique and symmetric angles with respect to it. Consequently, the cells take an elongated shape with the section of the nucleus that becomes elliptic with the long axis along the above directions as well. This fact well correlates with the observation that smooth muscle cells in the intima of arterial walls are oriented obliquely with respect to the vascular axial direction forming helical-like structures characterized by an angle with the longitudinal axis between 20 • and 40 • (Rhodin 1962;Shirinsky et al. 1989).
The reorientation dynamics in vitro is quite robust with respect to both cell type and experimental set-up. In fact, regarding the former aspect, fibroblasts, muscle-type cells, epithelial cells, endothelial cells, osteoblasts, melanocytes, mesenchymal stem cells, all respond in a similar way when periodically stretched. Regarding the latter aspect, the final result seems to be nearly independent from the applied frequency and amplitude, and from the mechanical characteristics of the substratum, with transitions when the corresponding values are smaller than some thresholds, i.e., too low frequencies, too small deformations, too soft substrata. On the other hand, the strain ratio in the two perpendicular directions turns out to be relevant, as well described by the experiments performed by Livne et al. (2014).
From the viewpoint of mathematical modelling, the first attempts to describe the phenomenon were based on a strain avoidance principle, consisting in the assumption that cells tend to reorient in the direction of minimal strain (Barron et al. 2007;Faust et al. 2011;Morioka et al. 2011;Wang 2000;Wang et al. 1995).
Successively, it was hypothesized that rather than minimal strain, the main reorientation direction tends to minimize stress (De et al. 2007(De et al. , 2008Livne et al. 2014). Therefore in these works, the evolution of the cell orientation θ is related to a linear elastic energy E through d dt In particular, Livne et al. (2014) modelled the ensemble of cells on the substratum as a linear elastic anisotropic material subject to biaxial strain and identified the equilibrium orientations θ eq formed by the cell major axis or of the stress fibers and the direction of stretching corresponding to minimal energy. In this way, they found a  Roshanzadeh et al. (2020). The top row refers to a case in which cells tend to orient at an oblique angle and the bottom row to the particular case in which ε yy = 0 (and then r = 0) for which cells tend to orient perpendicularly to the main streching direction linear relationship between cos 2 θ eq and a parameter quantifying the biaxiality of the deformation and the cell's anisotropic material coefficients. They also showed that in this parameter plane, data obtained using fibroblasts tend to align along a straight line. Starting from the observation that the experimental results holded true even for deformation ranges that make questionable the use of linear elasticity [they can go up to 30% (Faust et al. 2011;Livne et al. 2014)], Lucci and Preziosi (2021) proved that a generalization of the linear relationship found by Livne et al. (2014) also holds for a very large class of nonlinear constitutive orthotropic models. In the nonlinear framework, the squared cosine of the orientation angle is linearly dependent on a parameter which is the natural generalization of the one found by Livne et al. (2014), with a slope depending on a combination of elastic coefficients characterizing the nonlinear strain energy. A detailed bifurcation analysis is given. Also Lazopoulos and coworkers (Lazopoulos and Pirentis 2007;Lazopoulos and Stamenović 2006;Stamenović et al. 2009) employed a finite elasticity framework to describe stress fibers reorganization in strained cells, although they considered only uniaxial substratum stretching and addressed the problem using a non-convex energy, giving an explanation based on the co-existence of phases.
A viscoelastic model is instead proposed by Lucci et al. (2021) to explain why on the time scale of experiments the reorientation phenomenon does not occur for small frequencies, for instance, as a consequence of the reorganization of focal adhesions. A Maxwell-like force-deformation relation was also used by Chen et al. (2012) who focused on the dynamics of single stress fibers and on the reorganization of the attachment of focal adhesions to the substratum.
However, it must be noticed that for sake of simplicity most of the models mentioned above work in a deterministic framework, while, as in any biological process, randomness characterizes several aspects of the mentioned dynamics, such as the assembly and disassembly of stress fibers and of focal adhesions as well as the biochemical response inside the cell to such mechanical cues. Some of these aspects are considered in Hsu et al. (2009Hsu et al. ( , 2010, Kaunas et al. (2011) where the focus is on the stochastic evolution of radially oriented stress fibers around the nucleus when the cell is subject to static and cyclic stretch. De (2018) focused instead on the stochastic stretch-sensitive bond association and dissociation processes taking also into account the elasticity of the cell-substratum system to predict the orientation and stability of adhesion mechanisms.
From the experimental point of view, the visible result of such uncertainties is reflected in a spread in cell orientation, in the sense that the distribution of the orientations of the cells is not represented by a Dirac delta, but by smoother functions. Actually, the outcome of the experiments is naturally described using histograms and graphs reporting the distribution of the frequencies of cell orientations falling in a partition of angle ranges over [0 • , 90 • ] (see, for instance, Barron et al. 2007;Chen et al. 2018;Faust et al. 2011;Hayakawa et al. 2001Hayakawa et al. , 2000Livne et al. 2014;Mao et al. 2021;Neidlinger-Wilke et al. 2001, 2002Morioka et al. 2011;Wang et al. 1995;Wang and Grood 2000). The degree of spreading is not constant but depends on the amplitude and frequency of the imposed stretch. Specifically, it increases when decreasing amplitude and frequency.
The inclusion of some randomness allows the authors in Barron et al. (2007), Chen et al. (2018), Morioka et al. (2011), Wang et al. (1995 to compare the histograms obtained from the experiments with the curves obtained by the results of simulations of the orientation model that they propose. However, there, an analytical distribution function was not provided and the effect of stochasticity was not explored in detail. One of the first analitycal treatments of the problem of describing the probability density function of the orientations of the cells (its time evolution or, at least, the stationary state) is provided by Kemkemer et al. (2006Kemkemer et al. ( , 1999. They express the evolution of the orientation of a cell by means of an automatic controller, i.e. an ODE describing the temporal evolution of the single-cell orientation with an empirical forcing term that has the desired symmetry. They gain a stochastic differential equation (SDE) by adding diffusion, and obtain the evolution of the probability density function as a forward equation of the SDE. They can easily compute the stationary state of the resulting Fokker-Planck equation, represented by an exponential of a doubly-wrapped cosine, that is a Boltzmann-like distribution. In particular, they compare the analytical findings with experimental results and show that the Boltzmann-like distribution can describe cell orientations on curved substratums.
As a consequence, and as classically done in statistical mechanics when describing reorienting dipoles, many authors consider a Boltzman probability density function f kT , that is, as a matter of fact, coherent with the fact that the cells' orientation evolves according to (1). Then, all the efforts lie in the modelling of the energy E of the system and of its temperature T . For example, starting from their already mentioned works (De et al. 2007(De et al. , 2008, Safran and De (2009) describe the cell as a reorienting dipole subject to a periodic stretch and model the distribution of the orientations as a Boltzmann-distribution with a competition between the force determining the free energy of the dipole and the effective temperature. Faust et al. (2011) use this distribution assuming an E corresponding to the strain avoidance hypothesis. Also Mao et al. (2021) consider a Boltzmann-like distribution with an energy that is the sum of three contributions given by the work done by focal adhesions, by the pulling force, and by the elastic potential energy of bars in the tensegrity structure, that however presents a flaw. Driven by the aim of studying how peristalsic deformation affects the orientation of cells in the intestinal epithelial sheet, in a very recent work Gérémie et al. (2021), too, consider an SDE where the drift term is given by an elastic energy. They, then, determine the Fokker-Planck equation but they do not manage to retrieve a stationary probability density function, and they approximate it with a Gaussian distribution.
In the present work our aim is to determine the evolution and the stationary state of probability density functions describing the statistical distribution of the orientations of the cells subject to cyclic stretch. We shall do this using Fokker-Planck equations that we shall derive from microscopic stochastic processes taking into account the reorientation dynamics of cells in response to cyclic stretch. We shall then compare the probability density function with experimental results, in such a way that that the proposed microscopic process can be validated. As a first step, we shall consider a microscopic stochastic process ruled by the quadratic elastic energy proposed by Lucci and Preziosi (2021), that has already been investigated in a deterministic framework, leading to a good comparison with experimental results. As a second step, starting from the principle of minimization of the previous quadratic energy, we shall derive a microscopic stochastic process as a function of the actual rotation angles performed by the cell during the reorientation and caused by the cyclic stretch of the substratum. After deriving the probability density function describing the statistical distribution of the orientations of the cells that obey to this second stochastic process, we shall compare it with experimental results.
In order to do that, after recalling in Sect. 2 the mechanical background proposed by Lucci and Preziosi (2021), as a first step we shall model the evolution of the cell orientation by means of a stochastic differential equation in which the evolution of the direction is related to a general elastic energy plus a stochastic fluctuation (Sect. 3). In the same section the evolution of the probability density function is, then, classically obtained by means of a forward equation, namely a Fokker-Planck equation. We will find the stationary state and prove that it is an asymptotic equilibrium (Sect. 3.1). We will then show in Sect. 3.2 that using the elastic energy proposed by Lucci and Preziosi (2021) the results of the integration of the Fokker-Planck equation and its stationary state compare very well with the experimental results reported by Faust et al. (2011), Hayakawa et al. (2001, Jungbauer et al. (2008), Livne et al. (2014), Mao et al. (2021). In Sect. 4, we shall describe the process of reorientation as a discrete in time stochastic process that happens with a certain frequency. After showing that the same Fokker-Planck equation used in Sect. 3 can be obtained by performing a quasi-invariant limit of the Boltzmann kinetic equation describing the evolution of the statistical distribution of the orientations of cells, in Sect. 4.2 we will propose a different viewpoint that consists in modelling reorientation as a result of an internal optimal control problem activated by the cell. Finally, we compare the results of the integration of the derived Fokker-Planck equation and of its stationary state, obtaining an even better agreement with respect to the one obtained in Sect. 3.2.

Mechanical Backgrounds
We consider isolated cells that are seeded on the surface of a thin elastic substratum that is stretched biaxially. We define the x-axis along the direction subject to the maximum stretch. For sake of simplicity, we assume that cells behave elastically, are much softer than the substratum and strongly adhere to it, in such a way that the strain in the specimen is perfectly transferred to cells and is homogeneous. This translates in the fact that the strain tensor in the plane writes as E = diag(ε x x , ε yy ) = diag(ε, −r ε) where r is called biaxiality ratio. As sketched in Fig. 1, when stretched, cells internally develop stress fibers that link to the substratum via focal adhesions. The fact that these stress fibers tend to form along a certain angle with respect to the stretch direction, confers anisotropic characteristics to the system. Neglecting substratum deformability by the traction forces exerted by the cells, of adhesion remodelling, and of viscoelastic effect in cell behaviour [that are however considered in a deterministic fashion in Lucci et al. (2021), Xu et al. (2016)], we will describe the system through a general orthotropic elastic energy denoted by U that will be affected by the cell orientations.
Referring to Fig. 1, we will denote by θ the cell orientation angle with respect to the x-axis. Now, at variance with what happens during migration when the moving cell polarizes forming a head and a tail, in this case the internal structures of the cell aligned along θ and along θ + π are geometrically indistinguishable [see, for instance, the work by Roshanzadeh et al. (2020), Tondon and Kaunas (2014), Wang et al. (1995)]. As a consequence, these angles are also equivalent from the energetic point of view, i.e. one must have U(θ + π) = U(θ ). In addition, also the orientation of the axes is equivalent in the sense that it is observed experimentally that configurations θ and π − θ are equiprobable, as showed in Fig. 1. As a consequence, U(π − θ) = U(θ ). Therefore, U(θ ) is an even π -periodic function and we can work under the following symmetry requirements Most of the discussion that will follow is independent of the particular form of energy that is chosen provided that it possesses the symmetry properties above. However, using continuum mechanics arguments, it can be proved (see, for instance, the work by Ogden 2003) that an orthotropic elastic energy for a planar system with the symmetry properties in U 1 depends on θ through the square of its cosine and is characterized by material coefficients describing the response to stretching along the orientation direction (K ) and perpendicular to it (K ⊥ ) and the response to shear (K s ), in addition to possible mixed terms.
To be specific in the following we will neglect mixed terms and use the following form of elastic energy We notice that Eq.
We remark that in experimental works, observations of the orientations are reported over [0, π/2) when the parameters of the experimental setting correspond to cases (a) and (b), while the observations are reported over [0, π) when the parameters of the experimental setting correspond to cases (c) and (d). In particular, in the following we shall be interested in experimental settings that lead to scenarios (a) and (c) (Case 1 and 3). We want to highlight the fact that when dealing with scenario (a), experimentalists tend to represent data over [0, π/2) as, because of the aforementioned symmetry around the axis, the measure of θ ∈ [π/2, π) is reported in the histogram bin corresponding to π − θ ∈ [0, π/2). Working in a deterministic framework, on the basis of Lagrangian mechanics arguments, we can relate the evolution in time of the orientation angle θ with the changes in the virtual work L done by the stress acting on the cell due to stress fiber alignment. Considering an overdamped regime, which corresponds to neglecting inertial effects, we can then write where η is a viscous-like coefficient measuring cell resistance to internal rearrangement of stress fibers. In the elastic case Eq. (6) reduces to or where λ θ = η/K and we have put in evidence that the strain might be time-dependent. Referring to the work by Lucci et al. (2021) for a more detailed discussion, we here observe that the same equation is obtained for a viscoelastic Maxwell-like model in the limit of high frequencies ω with respect to the inverse of the viscoelastic relaxation time λ, i.e., λω 1. On the contrary, in the limit λω 1 viscous effects dominate and a term λω appears at the numerator (related to the appearance of a strain rate, i.e., ε(t)ε(t) instead of ε 2 (t)), so that the effective λ θ becomes λ θ λω . Considering that λ is of the order of one minute for both stress fiber and focal adhesion remodelling (Chen et al. 2013;Pasapera et al. 2010), one has that the transition from low to high frequencies occurs for ω about 0.01-0.1 Hz. In the work by Lucci and Preziosi (2021), Lucci et al. (2021) the authors perform simulations of the deterministic process (8), showing a very good agreement with experimental results that report the average orientation of cells subject to cyclic stretch.
At variance with the previous deterministic description and, as any biological process, cell reorientation is strongly affected by their stochastic behaviours. From the experimental point of view, then, this leads to a representation of the orientation state of the ensemble of cells in terms of mean, variance, and, whenever possible, frequency histograms, as discussed in the following (see Figs. 5, 6 and 7). In parallel, from the theoretical point of view, this leads to the need of determining a probability density function describing the statistical distribution of the orientations. For this reason, in the following, we will introduce a statistical mechanics approach.

Statistical Description of the Orientations of Cells under Bi-axial Stretch
In order to describe analytically the statistical distribution of the orientations of the cell, we introduce the probability density function f = f (t, θ) such that f (t, θ)dθ is the fraction of cells having orientation in [θ, θ + dθ ] at time t. As discussed above, the fact that cells have no identifiable head and tail, implies that if a cell is rotated by π , it is not possible to perceive a difference in cell orientation. Hence the angles θ and θ + π identify the same orientation. Therefore we shall deal with π -periodic probability density functions f , so that f (t, θ) = f (t, θ + kπ) ∀k ∈ Z. In addition, as a probability density function, f must satisfy Moreover, due to the symmetry related to the choice of the direction of the axes along the principal strain directions, the following property also holds F3: where F3 is also implied by the periodic character of the distribution function f .
With the aim of taking randomness into account, we may add a stochastic fluctuation to (7), where ξ is a Gaussian random variable with zero mean and unitary variance and σ takes into account the stochastic fluctuations due to uncertainties. The latter may then be more properly rewritten as an Itô process where dW t = √ tξ , being, then, W t a Wiener process. The Fokker-Planck equation describing the forward evolution of the probability density distribution f of the orientations of the cells that follow the dynamics (10) is then (Risken 1996) We observe that though in most experiments ε(t) = ε(1 − cos ωt), since we are interested in modelling the process of cell reorientation, as it is classically done in previously discussed elastic models, we will consider the mean strain ε over an oscillation period.
Introducing the nondimensional timet = tε 2 λ θ , then the Fokker-Planck equation whereσ 2 = σ 2 2ε 2 . This already puts in evidence that increasing the stretch amplitude decreases the dimensionless diffusion coefficientσ leading to a more focused response and more peaked distribution functions, and vice versa.
As already recalled, the inclusion of viscoelastic effects leads to the same results in the high frequency regime. On the other hand, in the low frequency regime, the dimensional analysis is modified because ε 2 is formally replaced by ε 2 λω. As a consequence, the effective dimensionless diffusion coefficient isσ 2 = σ 2 2λωε 2 , showing that when the imposed frequency decreases,σ increases leading to broader probability density functions.
We remark that Eq. (12) is similar to the one analyzed by Bastardis et al. (2008), and by Coffey et al. (2009) where the authors study a Fokker-Planck equation with a periodic potential that rules the rotational motion of a Brownian particle with inertial effects and that has the same symmetry properties as the elastic energy (2). In particular they extend Kramer's escape theory (Bastardis et al. (2008)) and treat a similar problem in the context of superparamagnetic relaxation of magnetic nanoparticles in 3D in Coffey et al. (2009).

The Stationary Equilibrium
Dropping the bars over f and t here and henceforth, if we denote by then the π -periodicity ofŪ and f , implies that In particular, thanks to the differentiability ofŪ, the stationary solution f ∞ of (12), coupled with an initial condition f 0 , satisfying F1, F2, F3 is found by imposing where the r.h.s. side is zero because of the boundary conditions (14). Thus, the stationary state of (12) is where C is a normalization constant. We observe that the maxima (resp. minima) of f ∞ (θ ) correspond to minima (resp. maxima) ofŪ. In particular, recalling that f is defined in [0, π), in Cases 3 and 4 there is only a maximum respectively in π 2 and 0. Therefore, in the former case, due to symmetry, the mean corresponds to the mode. A similar property can be obtained in the latter case working in the more convenient periodicity interval − π 2 , π 2 , otherwise the mean is trivially and misleadingly equal to π 2 . On the other hand, in Cases 1 and 2, f ∞ (θ ) is a bi-modal distribution in [0, π) with modes θ 1 eq , π − θ 1 eq and 0, π 2 , respectively. Actually, for the already mentioned symmetry reasons, usually, the range of angles used to report experimental data is the first quadrant [0, π/2), rather than [0, π) or [0, 2π). In this case, then, the notion of mean looses its informative role, especially with respect to the mode that, restricted to [0, π/2) is θ 1 eq in Case 1.

Remark 1
We observe that, if σ = 0, i.e. there is no stochastic fluctuation in (10), then the stationary state given by imposing (15) is a Dirac delta or a weighted sum of Dirac deltas centered in the stable equilibria orientations.
As usually done for the standard Fokker-Planck equation (Furioli et al. 2017), convergence to the stationary state can be studied by analyzing the monotonicity in time of various Lyapunov functionals of the solution. The typical one is the relative Shannon entropy where f , f ∞ : I ⊂ R → R + aree two probability densities. As periodic boundary conditions (14) hold, it is straightforward to prove (see Furioli et al. 2017) that the Shannon entropy monotonically decreases in time towards the stationary state, i.e.
Therefore, f ∞ is an asymptotic global equilibrium state.

Statistical Description and Comparison with Experiments
Usually, dealing with angles requires circular statistics and the definition of trigonometric moments (Mardia and Jupp 1999), e.g. the circular mean However, the symmetry properties of f in Case a) and Case b), that prescribe a bi-modal probability density function, would always lead to α = 0 and, therefore, θ(t) = π 2 . For this reason, we will use the following definition restricted to the first quadrantθ even because it better correlates with the definition of averagē that has been used in most experimental papers, where the 2 accounts for renormalization over [0, π/2). We will also use the coherent definition of variancē An index ∞ will identify the quantities above computed for the equilibrium distribution f ∞ . However, some remarks are needed. First of all, we observe that in general the average and the mode do not coincide, i.e.θ ∞ c ,θ ∞ = θ 1 eq . They obviously do when σ → 0. However, we will see numerically (see Fig. 4) that in most casesθ ∞ c =θ ∞ . In order to clarify this point, in Fig. 3 we plot the equilibrium distribution (16) over the interval [0, π) for different values of the parameters r andK s , α being fixed to the value α L = 0.794 determined fitting the data of the experiments by Livne et al. (2014). Then we varyK s and, from (3), set The positivity ofK ⊥ prescribes the compatibility conditioñ Therefore, we remark here that for all the figures presented most of the parameters are imposed by the experimental setting (r , α = α L ,K ,K ⊥ ,K s ). The values ofK are not in general given by experimentalists, butK = 1,K ⊥ is given by (21), and the only free parameters areK s , that must obey constraint (22), and σ . In particular we have observed thatK s does not influence the average of the distribution, as well as its shape, and we shall always considerK s = 0.7, while the greatest role is played by σ . We remark that in Fig. 3 and in all the other Figures, we have preferred to describe angles in degrees rather than in radians for a better readability and a more direct comparison with the statistical descriptions of the experimental results. In addition to the obvious observation that the diffusion parameterσ influences the spread of the orientations, other two facts linked to the presence of the diffusion stochastic term emerge explicitly and are put in evidence in Fig. 4: • unless for the symmetric case θ 1 eq = π 4 that is always obtained for r = 1 (see Eq.(4)), the average of the probability density distribution computed over 0, π 2 does not correspond to θ 1 eq , that is identified by the mode in the first quadrant, i.e. the maximum of the distribution function;  (4) and the mean θ ∞ over 0, π 2 (denoted by a circle in the figures), computed using (19) with f ∞ defined by (16) • the average of the probability density depends on σ and tends to the mode θ 1 eq (marked by ) when σ → 0 and to π/4, corresponding to a uniform distribution, when σ → +∞.
In Fig. 4 we also observe that the linear and the circular average at the stationary state coincide. Therefore, as experiments always consider the linear average, then in the following we shall make reference toθ . We remark that, for values of the average that are close to θ eq , there may be two different values of σ and therefore two different probability density functions that allow to recover the same averageθ (see Fig. 4). Therefore, at each time we shall determine the one that better reproduces experimental

Angle °F
ig. 5 Comparison of the evolution of the probability density function obtained by performing a Monte Carlo simulation of (10) with the experimental data reported in Hayakawa et al. (2001). In particular, ε = 20% and r = 0.4. Solution for θ 1 eq ≈ 61 • , σ ≈ 0.04, and λ θ ≈ 0.18 s. On the left, the thick bars in blue, red and yellow refer to the simulation results at t = 0, 1, 3 hours, respectively, while the corresponding lighter and thinner bars correspond to experimental datas. On the right, evolution of the probability density function.
results. It is evident that in Case 3 when θ 1 eq = θ 2 eq = π 2 , then it is more appropriate to use θ , rather thanθ .
With the aim of comparing the probability density functions with experimental results, we now focus on some papers reporting histograms of the percentage of cells in intervals of orientation angles. As in most cases esperimental data are given for θ ∈ 0, π 2 , we will restrict to the first quadrant. In Fig. 5 we compare the temporal evolution of the probability density distribution obtained by integrating (10) with a Monte Carlo approach with the experimental data reported in the work by Hayakawa et al. (2001) (the represented histograms) for ε = 20%, r = 0.4 and ω = 1 Hz, that implies that we are in a high frequency regime. In these experiments it is found that at t = 1 h the average orientation is 52.8 • , while at t = 3 hours, when more than the 80% of the cells are oriented at angles of 50 • -80 • , the average orientation is 62.02 • . Using (4) and α = α L the minimum of the elastic energy is obtained at θ 1 eq ≈ 61 • . In particular, we have run a Monte Carlo simulation of (10) with N = 10 6 cells and dt = 0.06 s, the initial probability density function is the uniform distribution over [0, π). We have recovered the probability density function as an histogram of the orientations of the simulated particles that, thus, approximates the solution to (12). In particular, we remark that the simulation is run over [0, π). We then restrict and renormalize f over [0, π/2) for comparison purposes with the reported histograms. Then, we calibrated σ in order to obtain a stationary state with average 62.2 • and that is closer to the histograms presented in Fig. 5 (left panel) and λ θ to replicate the time evolution of data. In particular, we set σ ≈ 0.04 that is such that θ = 62.2 • and the probability density function has the same height as the histogram at t = 3 h and λ θ ≈ 0.18 s. After 1 hour we have that the average orientation is 54.6 • and after 3 hours the average orientation is 62.04 • and 85% of the cells is oriented at angles of 50 • -80 • . In Fig. 5 (left panel) we plot both the histograms with classes' width of 10 degrees and, in the righ panel, the time evolution of the recovered probability density functions (that are histograms with classes width of 0.01 degrees). Focusing on the stationary distributions, Mao et al. (2021) report some experimental data in histograms over [0 • , 180 • ), changing the stretching amplitude (ε = 2%, 5%, 10%) and frequency (ω = 1 Hz, 0.001 Hz). In particular, they show that increasing values of both amplitude and frequency lead to more peaked distributions. In their case, r = 0 and the equilibrium orientation is perpendicular to the main stretching direction, i.e. θ 1 eq = 90 • . Trivially, due to symmetry, in this case the mode and the mean computed in [0 • , 180 • ) coincide, with σ, ε and ω determining only the variance of the probability density function. In Fig. 6 in order to replicate the data reported by the histograms by Mao et al. 2021, we plot (16) where we set the same σ = 0.2 and vary ε and ω. When ω = 1 (top row of Fig. 6), that corresponds to a high frequency regime, increasing the strain amplitude, coherently with the fact that σ 2 = σ 2 2ε 2 (so, it goes like ε −2 ) we have more peaked distributions that fit quite well the experimental distributions.
For ω = 0.001Hz since λω corresponds to a low frequency regime (it is λω = 0.1 if we take λ = 100 s), we useσ 2 = σ 2 2ε 2 λω . Also in this case, the distributions peak up when increasing the strain amplitudes and, the theoretical results compare well with the experimental results, in spite of the fact that we are not really using a viscoelastic model but only taking into account of viscoelastic effects through a modification ofσ that is valid in the low frequency regime. Comparing the results obtained for a fixed ε at the different ω's (for instance, the last column in Fig. 6) simulations give more peaked distributions for higher frequencies. Faust et al. (2011) report the results of some experiment characterized by an evaluated biaxiality ratio of r = 0.15. Assuming that α = α L , as also suggested in Livne et al. (2014), the minimum elastic energy and therefore the mode is obtained at θ 1 eq ≈ 79 • . They perform the experiment applying different stretching amplitudes, namely 4.9% (denoted as Case a 1 ), 8.4% (Case a 2 ), 11.8% (Case a 3 ), and 14% (Case a 4 ). We recall that in this case, at variance with the (symmetric) one in Mao et al. (2021), the mean changes with the strain amplitude that influencesσ (see second row in the table in Fig. 7). The means of the stationary distribution obtained by the simulation reported in the fourth row in the table closely follow the experimental ones. A slight difference is found for the standard deviation, expecially for larger amplitudes. Therefore, in Fig. 7 we compare their experimental results with the stationary probability density functions defined by (16) having average and standard deviation as computed from the histograms reported in Faust et al. (2011). In particular, we renormalize (16) over [0, π/2) for comparison purposes with the histograms that are reported in the work by Faust et al. (2011). We remark again that, as the average is close to θ eq in the presented cases, there may be two different values of σ and therefore two different probability density functions that allow to recover the same averageθ (see Fig. 4). Here, we have chosen the one that allows to better reproduce the histograms.

Kinetic Description
With the aim to get closer to the intrinsic dynamics followed by the single cell, in this section we will apply some classical tools of kinetic theory that, starting from the definition of the microscopic dynamics performed by the cells in the reorientation process, allow to derive the related mesoscopic evolution equation, such as (11). After going through the general procedure, we will then apply it to different microscopic rules. Then, in Sect. 4.2 we will discuss a different intrinsic dynamics that is probably performed by the cell, that through an optimal control argument drives them towards the most convenient orientation.  (2) in the cases a 1 , a 2 , a 3 , a 4 reported in the work by Faust et al. (2011) with applied strains listed in the table. In all figures we have r = 0.15 andK s = 0.7 that allowed to best reproduce the averages of the histogramsθ hist by varying σ in (16). The red circles represent the average circular orientationθ ∞ computed using (18). The black diamond represents θ 1 eq . We also computed the standard deviation of the histogramsd hist and the standard deviation v ∞ of the stationary state using (20).

Derivation of Kinetic Models from Discrete Random Processes
As a first step we formalize a microscopic discrete-in-time random process for describing the reorientation of cells. Let t ∈ [0, π) denote a random variable describing the orientation of a representative cell at time t. As typically done in kinetic theory (Pareschi and Toscani 2013), over a finite time interval t, we assume that a cell can change its main axis according to whether a reorientation occurs or not. We then express this discrete-in-time random process as where t is the random variable in [0, π) describing the new direction after a reorientation given the previous direction t , while T λ θ is a Bernoulli random variable which we assume to be independent of all the other variables appearing in (23), discriminating whether the direction changes (T λ θ = 1) or not (T λ θ =0) during the time interval t. In particular we set where naturally the necessary condition for T λ θ to be a random variable is t ≤ λ θ . Thus, the larger the time interval is, the higher the probability of having a reorientation is. The quantity t models the change of direction (if it happens) and it may be generally expressed as i.e., the new direction t is a function h λ,K of the previous orientation t and of the deformation parameters λ x , λ y , K , K ⊥ , K s , accounted for by the index λ, K . We shall assume h λ,K to be a regular function of its arguments, i.e. h λ,K ∈ C 1 ([0, π)), ξ is a standard gaussian random variable, i.e. ξ ∼ N (0, 1) satisfying ξ = 0, ξ 2 = 1, while the term mod(π ) models the fact that t is π -periodic. The aggregate description of the orientations of the cells can be obtained by determining the evolution of an observable quantity ϕ = ϕ(θ) defined on the phase space [0, π). Taking into account the rules (23) together with the assumed independence of T λ θ it is possible to see that the evolution of the probability density function f (t, θ) is (see "Appendix A" for a formal derivation) where, coherently with (25), θ is given by Equation (26) is a Boltzmann-type integro-differential equation. Choosing ϕ(θ) = 1 we readily obtain which means that the total mass of the agents is conserved in time by the interactions (27). Classically, the evolution of the statistical moments of f are obtained choosing ϕ(θ) = θ n , n = 0, 1, or circular moments may be recovered by setting ϕ(θ) = cos(θ ), sin(θ ).
As shown in "Appendix B", by an asymptotic procedure called quasi-invariant limit (see, for instance, Cordier et al. 2005;Furioli et al. 2017;Toscani 2006) based on a rescaled microscopic rule where γ 1, we can obtain in the limit γ → 0 a Fokker-Planck equation for the evolution of f In particular, if we want to model the new orientation of a cell that tries to minimize a potential energy U after a time interval dt we may observe that the discrete in time random process describing the evolution of the orientation t happens with frequency 1/λ θ and may be expressed by discretizing (10) over dt (where we consider the high frequency regime) and setting dt = γ Using the quasi-invariant limit procedure, we have the Fokker-Planck equation which is, as expected, the same as (11).

Reorientation as an Optimal Control Problem
In this section we want to introduce a new point of view consisting in modelling reorientation as a result of an internal control actuated by the cell that tries to minimize the elastic energy U. From the mathematical point of view, this approach consists in expressing reorientation rules like (27) starting from a control problem, in the sense that we assume that the cell changes its orientation by a rotation angle νψ opt where ψ opt is the angle that minimizes a certain cost functional J . At the kinetic level, this has been widely treated in recent literature, e.g. by Preziosi et al. (2021), Tosin and Zanella (2019), Albi et al. (2020), Dimarco et al. (2022), , . Therefore, we write where J is an energy functional defined as where the first contribution is a kinetic energy related to the control process, being ν a penalization coefficient, and the function g will be specialized later on. In order to determine the optimal control at each reorientation, we need to introduce a Lagrangian where χ ∈ R is the Lagrange multiplier associated with the constraint (30). The optimality conditions are then identified by the solution of Therefore, eliminating the Lagrange multiplier, the optimal value is implicitly identified by If we choose g = ε 2Ū , then Eq. (33) specializes to that, in general, allows to determine the optimal control only implicitly. In any case the reorientation rule (30) specializes into that in the limit of small ν used for the grazing limit and adding the stochastic term is equivalent to (29) and leads again to (12). In order to explicitly determine the control, we can, instead, more classically take a quadratic form for g where, assuming to work in Case 1, with p(θ ) a non negative and continuous function defined on [0, π) that satisfies p(θ 1 eq ) = 1 p(π − θ 1 eq ) = 0, p(0) = p(π/2) = p(π ) = 1/2, p (θ 1 eq ) = p (π − θ 1 eq ) = 0 in such a way thatθ (θ 1 eq ) = θ 1 eq andθ(θ 2 eq ) = θ 2 eq .
Therefore, the choice of g given by (34)-(37) has the same extrema θ 1 eq and θ 2 eq as U . The latter models, essentially, the fact that if θ is already close to an equilibrium orientation, then it is more likely not to change. In particular, we shall consider a second order polynomial satisfying the previous conditions. In this case one can explicitly solve (33) and determine ψ opt = − ε 2 1 + νε 2 (θ −θ), and therefore the reorientation rule (30) becomes Adding a stochastic fluctuation weighted by σ c we have This rule implies the fact that at each reorientation the cell will activate a control to reach a better orientation that is given by a rotation of γ ε 2 (θ(θ) − θ) (plus a white noise). This process will stop when the cell has oriented along the stable equilibria, because of the choice (34)-(37). In the symmetry points θ = 0, π/2, π the cell has the same probability (= 1/2) of reorienting either towards θ 1 eq or θ 2 eq = π − θ 1 eq . As illustrated in Sect. 4.1, in this case the quasi-invariant direction limit procedure leads to the following Fokker-Planck equation that can be coupled with boundary conditions F3. Therefore, the stationary state is given by where C is the normalization constant. This distribution has actually mode θ 1 eq and π − θ 1 eq in [0, π), thanks to the choice (34)-(37) and average depending on the value of σ c .
In Fig. 8 we compare the stationary distribution (41) with the experimental data by Faust et al. (2011), as in Fig. 7. Setting σ c in such a way thatθ l of (41) is the same as in  Fig. 7, in all figures we have r = 0.15 andK s = 0.7 that allowed to best reproduce the averages of the histogramsθ hist l by varying σ c in (41). The red circles represent the average circular orientationθ ∞ l computed using (18) with (41). The black diamonds represent θ 1 eq . We also computed the standard deviation of the histogramsd hist and the standard deviation v ∞ of the stationary state using (20) with (41). We also superpose (16) withŪ given by (2) as reported in Fig. 7 the work by Faust et al. (2011), we find that the microscopic rule (39) allows to recover probability density functions (41) that are better than those in Fig. 7. The prediction of the standard deviation, reported by the fourth line of the table in the two figures, shows that those of (41) are slightly closer to the linear standard deviation reported by Faust et al. (2011). We remark that the values of σ and σ c are very different, and this is due to the fact that the rule (10) expresses the variation of θ in terms of its derivative and of the elastic energy, while (39) expresses the variation through a rotation angle that the cell performs during a reorientation. Fig. 9 Temporal evolution of the mean of the orientation distribution. In (a) ω = 1.2 Hz and ε = 10% as reported by Livne et al. (2014). In addition, λ θ = 6.6 s and σ c = 0.7. In (b) and (c) ω = 2 Hz and ε = 8% as reported in Jungbauer et al. (2008). In addition, λ θ = 6.6 s and σ c = 1.6. After 3000 s stretching stops and cells tend to reorient uniformly. The standard deviation (one confidence interval) of the angle is also given in (b). In c the same mean is reported in terms of its cos 2θ for a more direct comparison with the work by Jungbauer et al. 2008. Green squares correspond to the experimental results reported by Jungbauer et al. (2008).
Focusing on the temporal evolution of (40) in Fig. 9(a) we report the results obtained performing a Monte Carlo simulation of (23), (24), (39) with N = 10 6 particles, γ = 10 −2 , as done in a different context for example by Loy and Tosin (2021). In fact, equation (40) is derived as the quasi-invariant limit of a Boltzmann like equation (see "Appendix B") with microscopic rule (39), that is derived from (23)-(24)-(25) in the limit of large N and small t. In particular, we choose the data of the experimental results reported by Livne et al. (2014) where ε = 10%, λ θ = 6.6 s and ω = 1.2 Hz, corresponding to a high frequency regime, and we set σ c = 0.7 so that the average orientationθ l of (41) withσ 2 c = σ 2 c 2ε 2 is the same as reported in Livne et al. (2014). The qualitative behaviour corresponds to that reported in Livne et al. (2014). In particular we find that the rotation time is λ θ /ε 2 as stated in Livne et al. (2014).
Eventually, we want to replicate the experiment proposed by Jungbauer et al. (2008), where the authors stop stretching at a certain time and record the recovery phase towards a uniform distribution. To this aim, in Fig. 9b, c, we simulate (39) with N = 10 6 elements, γ = 10 −2 and using the parameters corresponding to the experiment reported in the work by Jungbauer et al. (2008): the stretch is imposed only for 3000 s, while ε = 8% and r = 0.194. After 3000 s, ε = 0. We choose the same reorientation time as found in the work by Livne et al. (2014), i.e. λ θ = 6.6 s, for the whole dynamics. Also in this case the behaviour corresponds to that reported in the work by Jungbauer et al. (2008) (green squares corresponds to the experimental results reported in Jungbauer et al. (2008)).

Discussion
In order to describe the dynamics of cell reorientation under stretch, we have proposed a class of Fokker-Planck models for the evolution of the statistical distribution, i.e. the probability density function, paying particular attention to their link with the microscopic models. In particular, we have considered a stochastic microscopic process (10) in which the cell tends to minimize the elastic energy U. The model is able to describe both the evolution and the stationary state of the probability density function over the orientations of the cells, which can be determined explicitly from as the stationary state of the Fokker-Planck equation relative to the SDE (10). The results compare well with several independent experiments (Faust et al. 2011;Hayakawa et al. 2001;Mao et al. 2021) showing the flexibility of the model.
In Sect. 4, we have used a well known procedure that allows to recover Fokker-Planck equations from microscopic stochastic discrete in time processes, through classical tools of kinetic theory. We have shown that by means of this approach it is possible to recover the Fokker-Planck equation (12) thanks to an appropriate choice of the microscopic rule for the evolution of the orientation angle. Then, thanks to the optimal control problem we have obtained a rule that is expressed as a function of the rotation angle performed by a cell during a reorientation. Also in this case the results compare well with several independent experiments (Faust et al. 2011;Jungbauer et al. 2008;Livne et al. 2014).
At present, the microscopic dynamics determining the drift term in the Fokker-Planck equation is defined according to biophysically sound qualitative arguments. In the future, the close link between the microscopic and the mesoscopic model shown here can be exploited, on the one hand, to better calibrate the model with respect to experimental data, and, on the other hand, to describe the microscopic mechanisms starting from measurements on the behaviour of single cells, whenever these data will be experimentally available. Moreover, the advantage of the microscopic rule (39) is that it is expressed in terms of rotation angles and is thus more amenable to possible extensions in order to include superposing effects that can be considered when dealing with cells seeded on a substratum, for example on collagen, that are subject to cyclic stretch, such as contact guidance and steric hindrance (Ristori et al. 2018). Moreover, the present framework may be extended to describe a three dimensional environment, by considering a second angle and its microscopic dynamic and a probability density function that depends on the two rotation angles.
we get which, by integration by parts, and recalling the compactness of the support of ϕ, can be recognised as a weak form of the following Fokker-Planck equation