Modeling cells spreading, motility, and receptors dynamics: a general framework

The response of cells during spreading and motility is dictated by several multi-physics events, which are triggered by extracellular cues and occur at different time-scales. For this sake, it is not completely appropriate to provide a cell with classical notions of the mechanics of materials, as for “rheology” or “mechanical response”. Rather, a cell is an alive system with constituents that show a reproducible response, as for the contractility for single stress fibers or for the mechanical response of a biopolymer actin network, but that reorganize in response to external cues in a non-exactly-predictable and reproducible way. Aware of such complexity, in this note we aim at formulating a multi-physics framework for modeling cells spreading and motility, accounting for the relocation of proteins on advecting lipid membranes. We study the mechanical response under compression/extension of an assembly composed of 8 helical rods, pin-jointed and arranged in pairs with opposite chirality. In compression we find that, whereas a single rod buckles (a), the rods of the assembly deform as stable helical shapes (b). We investigate the effect of different boundary conditions and elastic properties on the mechanical response, and find that the deformed geometries exhibit a common central region where rods remain circular helices. Our findings highlight the key role of mutual interactions in the ensemble response and shed some light on the reasons why tubular helical assemblies are so common and persistent.


Introduction
Receptors dynamic along cell membrane is a key factor in several biological phenomena, as for angiogenesis, tumor metastasis, endocytosis and exocytosis. Angiogenesis is a multistep process in which endothelial cells are affected by several extracellular stimuli, including growth factors, extracellular matrix, and parenchymal and stromal cells. In this process, growth factor receptors as well as adhesion receptors convey the extracellular signaling in a coordinate intracellular pathway promoting cell proliferation, migration, and their reorganization in active vessels [1]. Integrins are a family of cell adhesion receptors that support and modulate several cellular functions required for tumor metastasis. They can directly contribute to the control and progress of metastatic dissemination. During tumor development, changes in this family of receptors impact upon the ability of tumor cells to interact with their environment and enable metastatic cells to convert to a migratory and invasive phenotype. Integrins regulate each step of the metastasis and affect tumor cell survival and interaction with changing environments in transit from the primary tumor to distant target organs [2]. Receptormediated endocytosis is a process by which cells absorb metabolites, hormones, proteins and, in some cases, viruses by the inward budding of the plasma membrane (invagination). This process forms vesicles containing the absorbed substances and is strictly mediated by receptors on the surface of the cell [3].
Whereas uncountable papers have been published on the biology of cells spreading, motility and the relocation of proteins on advecting lipid membranes, the mathematical modeling definitely lags behind experiments and overall received much less attention. Although nowadays a widespread literature in mechanobiology exists [4], the relocation of proteins and their interaction with the reorganizing cytoskeleton in the biological phenomena mentioned above is still an ongoing research topic, let alone the formulation of efficient algorithms and computational solvers for three-dimensional simulations [5].
In this note, we attempt at defining a multi-physics scheme for the modeling of cells spreading, motility and the relo-cation of proteins on advecting lipid membranes, framing the mathematical setting within the mechanics and thermodynamics of continua [6], stemming from seminal works [7][8][9] and accounting for recent literature, either connected to the endocytosis of virus in human and animal cells [10][11][12] or ligand-receptor mediated raft formation [13], chemotaxis [14], surface-associated caveolae mechanotransduction [15]. The general framework illustrated in this note applies to growth and remodeling, too, falling within the category of theory of finite growth according to the terminology defined in [16].
The paper is designed as follows. After a nomenclature of the main symbols and the definition of operators in a Lagrangian setting, the paper focuses in Sect. 3 upon the relocation and reaction of receptors on a lipid membrane that advects. The topic is purposely presented in a broad sense, in order to be applicable to several possible receptorsligands interactions: specific applications-carried out in [17,18] and shortly illustrated here in Sect. 7-deals with the relocation of vascular endothelial growth factor receptors and integrins during endothelial cell adhesion and spreading. In spite of the generality, Sect. 3 is self-contained and includes the description of Reynold's theorem on a surface that advects, of the equations that rule proteins transport on an advecting lipid membrane, and eventually of the receptors-ligand interactions, in form of chemical reactions, that take place concurrently with relocation. A rather similar approach has been taken in Sect. 4, which concerns the relocation and reaction of actin to form biopolymers within the cytosol. The mechanical evolution of the cell is discussed afterwards in Sect. 5: besides stating the classical balance laws (of linear and angular momentum), the Sect is accompanied by an extensive discussion on boundary conditions, aimed at showing that Neumann type of conditions, due to electrostatic interactions, are most likely not responsible for cell spreading and motion in view of the modest amount of energy involved in those interactions compared to the bulk energy of a cell. We concluded therefore that spreading is a result of extensional and contractile forces exerted by pseudopodia and the cytoskeleton machinery [19]. Those forces have been investigated further in Sect. 6, where the thermodynamics of receptors motion on the membrane was studied at first up to the constitutive theory and the receptors-ligand interactions kinetics. The analysis of the thermo-chemo-mechanics of cells is the last Sect of this work: in it, we highlight the role of strain and stress decompositions in order to model cell adhesion, protrusion, and contractility. A bibliographic review is presented in a rather extensive paragraph, showing various approaches pursued in the literature to cover the multiscale scenario of cell viscoelasticity and identifying missing pieces within the theoretical framework that we set in the present note.

Definitions
Denote with Ω(t) a volume that advects, and with ∂Ω(t) its surface. A point x ∈ Ω(t) is defined as the image of a point X in a reference configuration Ω R through a smooth function χ(X, t) termed motion. Following [6], we will name deformation the snapshot of a motion at a fixed time t: The deformation is assumed to be a one-to-one map. In addition, denoting the deformation gradient with the requirement J = det [ F ] > 0 holds. Define on the surface a part P(t) ⊂ ∂Ω(t) as in Fig. 1, and consider a scalar function f (x, t) with x ∈ P(t). Denote with v adv (x, t) = dx/dt the velocity of advection at location x and time t; such a velocity has an arbitrary direction, i.e. it is not necessarily tangent to ∂Ω(t).
The Frenet-Serret reference frame at a generic point y ∈ ∂P(t) is defined as in Fig. 1, in terms of the two unit vectors t (y, t) (tangent) and t ⊥ (y, t) (normal). The vector n(y, t) (binormal) is here taken of non-unit length, being the imagine in Ω(t) of a unit vector n R in the reference configuration Ω R , by means of the contravariant transformation On the other hand, the following covariant transformations hold: with the obvious implication that t R and t ⊥ R are not unit vectors. The Frenet formulae holds, namely: where κ denotes the curvature and τ the torsion.
The projected gradient operator of a scalar field f on a surface P is defined as follows Note that x ∈ P(t) implies X ∈ P R . (b) Frenet frame at point y ∈ ∂P(t) and the normal vector n at point x ∈ P(t) in the current configuration, whereas in the reference configuration it reads The projected divergence operator of a vector field v, which has an arbitrary direction, on a surface P is defined as follows in the current and reference configurations, respectively. Tensor l is the gradient of v, l = ∇ [v]. Note that l in Eq. (2a) can be replaced by its symmetric part d = sym [ l ], since for any skew-symmetric tensor w it holds n · wn = 0. Alternative forms for the projected divergence operators are Provided sufficient smoothness, the divergence theorem holds also for advecting membranes, in the form: The proof of this theorem, as well as for all other theorems not explicitly stated in this paper, can be found in [20].

Reynold's theorem on a surface that advects
Reynold's theorem on P(t) reads as follows [20]: where v adv (x, t) is the velocity of advection at location x and time t. By taking f = 1, Eq. (5) depicts the area evolution of P(t)as It is intuitive that advection with velocity in the tangent plane has the potential of modifying the surface area, however even v adv (x, t) ∝ n(x, t) can do so, as for the homothetic expansion of a rubber balloon. Reynold's theorem (5) can be also restated as (6) and is a restriction on surfaces of the classical Reynold's transport relation on volumes ( see [6], Sect. 16 among others).

Mass balance in the current configuration for a convecting species
Consider a generic species a at a point x on the surface ∂Ω(t).
Species a convects with velocity v a (x, t). The latter entails the dragging, or advection, velocity v adv (x, t) and another velocity that is due to many possible physics, as for diffusion or migration. If internalization of species from the membrane is not allowed, the net velocity v a − v adv lays in the tangent plane of the membrane and Since species are modeled on a membrane, which is a twodimensional manifold, the surface density ρ a of species a measures the mass of the species per unit surface. The density flux vector of species a, denoted withh a , is the product of the surface density times the net velocity of species a, i.e.
Define on the surface a part P(t) ⊂ ∂Ω(t) as in Fig. 1. The flux of species a across the boundary ∂P(t) is and the mass balance of species a in the advecting configuration P(t) reads where s a (x, t) is the surface mass supply 1 of species a. By means of the divergence theorem (4) and of Reynold's transport theorem in the form (6), balance law (9) becomes Since it holds for all P(t), it eventually localizes as This formulation of the mass conservation law has been considered also in [21]. The mass balance can be finally written in terms of surface molarity c a (in moles or molecules per unit surface), by division by the molar or molecular mass (m a ) of species a. By denoting with c a = ρ a /m a , s a = s a /m a , and h a =h a /m a the local balance (11) becomes

Mass balance in the reference configuration for a convecting species
The mass balance (12) can be rephrased in the reference configuration at point X and time t. To this aim, define the reference molarity of species a as the reference flux vector h a R (X, t) and the reference mass supply s a R (X, t) as respectively, where [6,22]: The referential form of the mass balance (12) can be derived from the mass balance in the form (9), and reads For the sake of brevity, the proof has been here omitted, interested readers may find it in [20].

Relocation and reaction
The association and formation of a protein complex follow a two-steps mechanism; the formation of an encounter complex, in which previously free proteins show few specific interactions and assume many orientations, and the evolution of the encounter complex in the final complex. The bindingunbinding interaction accounts for both mechanisms [23]. Coefficients k f and k b are the kinetic constants of the forward and backward reactions, respectively. The rate of reaction (17), denoted with w (17) and measured in [ mol m 2 s ], quantifies the net formation of (C) on the advecting membrane as the difference between the forward and backward reactions. Equation (16) shall be extended to account for the reaction (17) and tailored to species a = R, L, C. Receptors (either free or bound into the complex) are distributed along the membrane together with other lipid species and proteins. They are assumed to freely move laterally, effects due to steric hindrance are not accounted for. The amount of proteins per unit area that can be placed at a membrane location x is thus limited by the actual size of the protein itself. This evidence ushers the definition of a saturation limit for the species, c max a (x, t). During their life, cells and their membranes undergo major macroscopic mechanical deformations. Studies on the red blood cell [24] suggest that the membrane deformation occur at constant area, but this evidence does not appear to be supported by experiments in endothelial cells during spreading [19]. Individual protein and phospholipid can easily move laterally within the membrane, which results in a very low shear stiffness. The fluid mosaic model [25] captures this evidence, adding a questionable high resistance to areal expansion. Indeed the mechanisms that are in charge of areal expansion during cell spreading are complex and involve the microstructural topology 2 of the membrane (as for flattening of invaginated membrane domains [26], i.e. the role of the caveolae as membrane surface repository readily made available for fast geometrical evolution as during filopodia extension). The structure of the lipid membranes, however, induce to suppose that the saturation concentration c max a (x, t), i.e. the maximum number of moles or molecules per unit area for any species a, remains unchanged in time in the current config-uration. This choice in turn entails that the number of moles or molecules per unit area in the reference configuration is not constant and evolves in time following Eq. (13), i.e.
Accordingly, the value of the non-dimensional ratio between the concentration of species a and its amount c max a at saturation, (19) in the current configuration remains unchanged in the reference configuration The kinetics of reaction (17) is modeled as for ideal systems via the law of mass action [27] w (17) At chemical equilibrium, as w (17) = 0, the concentrations obey the relation (22) which defines the constant of equilibrium K (17) eq of reaction (17).
Far from the saturation limit, (1 − ϑ a ) ∼ 1 for all a. Accordingly, the mass action law (21) simplifies as w (17) once the new constants The diffusion of receptors and the viscous evolution of the cell during adhesion and migration appear to be much slower than the interaction kinetics, i.e. the time required to reach chemical equilibrium is orders of magnitude smaller than the time-scale of other processes. For this reason, thermodynamic equilibrium may be invoked in place of a transient evolution, thus inferring the constraint w (17) = 0 to the concentrations of species at all times. Far from saturation, equating (23) to zero implies that having denoted with α the following constant: In view of identity (24), the two concentrations c R and c L describe the problem in full, and the concentration of the complex can be deduced a posteriori.
In vivo experiments show that the complex molecules usually have a much smaller mobility than receptors, perhaps induced by their size. For in vitro experiments [17,18,28], ligands are prevented to flow onto the substrate: given that complex molecules result from the interaction with immobile ligands, they are macroscopically steady as well. Since receptors move along the membrane, reaction (17) traps mobile receptors and vice-versa [29]. In this work, analogously to [30], ligands and complex are assumed to be motionless, i.e.
The reaction rate w (17) (x, t), being a mass supply, shall transform as s a (x, t) according to Eq. (14). The invariance of ϑ a with the configuration and the analysis of the mass action law (21) imply that the forward and backward "constants", which encompass the dimensionality of w (17) (x, t), are not actually constants in the reference configuration. They rather change with time and with point X according to with j(X, t) as in (15). The equilibrium constant in the reference configuration, being the ratio of k f R and k b R remains independent upon the configuration. Eventually, the mass action law (21) in the reference configuration writes In view of all considerations made so far, the local form (16) of the mass balance specify as follows (omitting the dependency upon X and t ): where the receptors flow. The supply s R R accounts for internalization or generation of proteins: it is the amount of receptors that are generated within the cell and reach the membrane or that internalize. It can be related to the change in the membrane area through a parameter κ R R as At all points at which ligands and receptors do not interact, the reaction rate w (17) R vanishes. Equation (29b) is rather defined in the location where ligands stand. In vitro, a given amount of ligands (which can be thought of as the initial condition of Eq. (29b)) are spread upon a microscope slide. Finally, Eq. (29c) is defined in the contact zone between the cell and the slide where reaction (17) takes place.
It is convenient to rephrase Eq. (29b) in terms of the "ligands made available for the reaction" in place of the "ligands spread on the slide". The former ligands are the ones "felt" at a point on the membrane as the distance from such a point and the substrate, where ligands are spread out, becomes sufficiently small.
Such a distance can be understood as a cutoff, within which the formation of an encounter complex, C * , becomes possible as a consequence of diffusion, as made clear in [23,[31][32][33]. Despite the size of the cutoff distance remains inaccurately estimated, it was established to be on the order of tens nanometers [7,23]. It arises form the interplay of attractive and repulsive forces between either two cells or a cell and a substrate. Indeed, negative electrical charge carried by cells generates repulsive electrostatic forces-repulsive barrier -which is further enriched by an additional resistance provided by the compression of the glycocalyx proteins. Rather, electrodynamic van der Waals forces are expected to be attractive [23]. Both van der Waals and compressive forces are characterized as non-specific long ranged forces, whereas cell adhesion is generally mediated by the specific short ranged receptor-ligand interactions, which can cause cell adhesion much more tightly than the non-specific electrical forces [23,30]. Cells separated by a distance less than, or equal to, the cutoff distance should form a zone of adhesion with the substrate by means of local fluctuations in receptors density, so that small regions of increased density can penetrate through the resisting potential to react with the source of ligands on substrate [7].
This point of view, which corresponds to the picture of tight receptor-ligand bond as a set of weak non covalent physical interactions [34], is made explicit by a supply function s L R , that vanishes at long ranges and rapidly reaches the initial concentration of ligands available for the reaction at short distances The ligand supply s L R (X, t) becomes available for the reaction during the spreading of the cell. It seems to be logically related to: i) a gap function between the substrate rich in ligands and the cell membrane in the current configuration; ii) a lag in time, namely a point-wise function of an internal variable that activates when the gap function is below some threshold and is related to the chemical kinetics of the binding-unbinding reaction (17). In this form, all three equations (29a), (29c), (31) can be written on the membrane X ∈ ∂Ω R . Assuming that the time scale of the chemical reaction is much faster than other processes, the concentrations of species may be governed by thermodynamic equilibrium at all times. The concentration of complex c C R relates then to the others by the equation w (17) = 0, which leads to Eq. (24) in the current configuration. Making use of mapping (13), Eq. (24) relates the concentration of complex in the reference configuration c C R to the concentration of ligands and receptors in the same configuration c L R , c R R as follows with constant α defined in Eq. (25). Transformation (32a) is consistent with the assumption (18) made on how saturations transform.
In conclusion, exploiting identity (32a), the two concentrations c R R and c L R fully describe the problem in the assumption of infinitely fast kinetics, whereas the concentration of the complex can be deduced a posteriori. The two governing equations descend from Eqs. (29) and read: Equations (32), with associated initial conditions and Dirichlet-Neumann boundary conditions define the relocation of receptors that undergo binding-unbinding reactions on the reference configuration of a membrane that advects. These are balance equations and as such hold for any constitutive behavior for the mass flux. These equations are coupled to the mechanical evolution of the cell (i.e. adhesion, spreading, migration) through the function s L R (X, t), which "transfers" ligands on the membrane according to the geometry of the cell.

Relocation and reaction of actin to form biopolymers
The extensive mathematical description made in Sect. 3 will guide the modeling of the relocation and reaction of actin to form biopolymers in the cytosol, which will be summarized here in a shorter shape. Biopolymers are composed of actin, a protein termed globular or G-actin in its monomeric form and F-actin when it forms filamentous polymers-see Fig. 3(a). In turn, actin filaments can bundle to form stress fibers, or cross-link to form polymer networks that allow the movement of the cell-see Fig. 3(b). Polymerization is usually triggered by extracellular signals. In the case of cell locomotion, for instance, the cell extends finger-like protrusions by which the cell "feels" the surrounding surface. As done in [35], the precise details of the signaling pathways are here ignored. Rather, the level of signaling is assumed given in the reference configuration by a function that accounts for the location of discrete signaling points y i in the surroundings emitting signals of intensity γ i at time τ i ; θ is the decay constant of the signal. This approach in modeling the external stimulus is similar to the membrane activator in Ref. [36]. The transduction of the signal results in the polymerization of the actin filaments and their cross-linking or bundling. The formation of single actin filaments can be modeled as a bimolecular reaction, as in [37]; in this note, the biopolymer turn-over will be described at a larger scale, involving the interplay between fundamental units and stress-fibers or pseudopodia, in the form ( 3 4 ) with F denoting either one of the two biopolymers. The network or fiber formation rate of reaction (34), denoted with w (34) , is influenced by mechanical stresses: stress fibers stability is favored by tension, for instance. For this reason, the stress tensor enters the chemical potential and the dissociation reaction of biopolymers. The kinetics of reaction (34) is modeled via the law of mass action, properly extended to account for signaling: having already discussed the meaning of the ratio ϑ in Eq. (19). Function D accounts for the role of the stress in the dissociation of biopolymers, see for instance [35].

Mass transport in the cytosol
Consider a generic species a at a point x in the cytosol Ω(t).
The mass balance of species a in the advecting configuration localizes as withh a and v adv defined earlier in Sect. 3.2.1, ρ a is the density of species a. The mass balance can be restated in terms of molarity c a (in moles or molecules per unit volume), by division by the molar or molecular mass (m a ) of species a. By denoting with c a = ρ a /m a , s a = s a /m a , and h a =h a /m a the local balance (36) becomes The latter can be rephrased in the reference configuration at point X and time t. To this aim, define the reference molarity of species a as the reference flux vector h a R (X, t) and the reference mass supply s a R (X, t) as [6] h a R = J F −1 h a ( x(X, t ), t ), respectively. The reaction rate w (34) (x, t), being a mass supply, shall transform according to Eq. (39)b. The invariance of ϑ a with the configuration and the analysis of the mass action law (35) imply that the forward and backward "constants", which encompass the dimensionality of w (34) (x, t), are not actually constants in the reference configuration. They rather change with time and with point X according to The ratio k f R /k b R remains independent upon the configuration. The referential form of the mass balance equations eventually reads As for the complex molecules, filaments usually have a much smaller mobility than monomers and might be assumed to be motionless, i.e.
The diffusion of monomers appears to be much slower than the interaction kinetics and the concentrations of species may be governed by thermodynamic equilibrium at all times [38]. The concentration of filaments c F R relates then to the monomers by the equation w (34) = 0, mediated by the local amount of signaling and stress. Equations (41), with associated initial conditions and Dirichlet-Neumann boundary conditions define the relocation of monomers that undergo polymerization reactions in the reference configuration.

Mechanical evolution of the cell
Based upon the selection of the mechanisms that are supposed to govern the structural response of the cell, the balance laws of linear and angular momentum come out. Literature provides two basic approaches, whether the structural functions are demanded entirely to the cell membrane [39][40][41][42][43] or to the development of a cytoskeletal structure within the bulk of the cell [9,35,[44][45][46][47][48][49][50]. The influence of curvature on the elastic stiffness of the membrane appears to be related to the size of the cell [51] and seems to be negligible for endothelial cells of diameter ∼ 10 µm. These two evidences lead to consider the reorganization of the cytoskeleton through a network of actin and intermediate filaments and microtubules the main responsible for the mechanical response of endothelial cells [52], coupled to a passive behavior dictated by the viscosity of the cytosol as in [9,35,50]. Accordingly, balance of linear and angular momentum will be formulated for the bulk of the cell rather than the membrane. Forces in continuum mechanobiology are described spatially by contact forces between adjacent spatial regions (as for the forces exchanged by the substrate and the cell during adhesion), surface forces exerted on the boundary of the cell by the environment (as for the receptor-ligand attractive interaction [23,53] and repulsive electrostatic interactions), body forces exerted on the interior points by the environment (as for the gravity or pseudopodia forces that preside migration). Contact and surface forces, acting on ∂Ω(t) will be denoted henceforth with t(x, t) whereas body forces will be denoted with b(x, t). Their referential counterparts will inherit the subscript R .
Throughout the rest of the paper we will neglect inertia forces, although some authors [54] pinpointed the role of inertia forces during migration. Accordingly, the balance of linear and angular momentum, which are assumed to hold at each time for all spatial regions Q(t) ⊆ Ω(t), read: with r denoting the position vector with respect to an arbitrary pole. Classical arguments of continuum mechanics lead to localize Eqs. (43) in the reference configuration, in terms of the (first) Piola stress tensor P and of the body forces measured per unit volume in the reference body b R (X, t) = J (X, t) b(x(X, t), t).
The referential local form of the balance of linear momentum reads The first Piola stress tensor P must satisfy the local angular momentum balance Contact and surface forces are boundary conditions for problem (44a). They emanate from electrostatic long or short range interactions, from receptor-ligand adhesion forces, as well as from contact tractions after adhesion. A vast literature [37,55,56] has been devoted to quantify the forces involved in these interaction mechanisms. It emerges that uncertainties remain in the establishment of realistic values for attraction forces, not surprisingly due to the complexity of the required experimental tasks.

Thermodynamics
The quest of the right thermodynamic principles in mechanobiology is, on one hand, far from being understood and, from a wider perspective, it paves the way to boundless questions of philosophical and ethical nature, as for the establishment of a thermodynamics of life [57], which fall completely out of the scope of present paper. Major accomplishments have been recently achieved [58] in formulating fresh concepts that deviate from classical results of thermodynamics of non equilibrium. In this scientific area, which is nowadays flourishing, new fundamentals assertions are expected in the years to come. Being aware that classical formulations of non equilibrium thermodynamics [27] may not be able to capture some principles of mechanobiology that rule the dynamic of receptors-as for the homeostatic constraint [59], we are prone to deepen our formulation in future studies.

Thermodynamics of receptors motion on the membrane
Following [29], receptors motion on the lipid membrane is thermodynamically described by energy and entropy balances, imposing as usual that the internal production of entropy cannot be negative. After the definition of the referential specific Helmholtz free energy per unit volume ψ R , taken as a function of temperature and concentrations, ψ R T , c R R , c L R , c C R the entropy imbalance in the Clausius-Duhem form is derived. Standard arguments-the so called Coleman Noll procedure-finally allow identifying the following thermodynamic restrictions for chemical potentials μ and entropy η R . Assuming further that relocation of receptors take place in thermal equilibrium conditions, the so called Clausius-Plank inequalities apply: R the chemical affinity of reaction (17). A strategy to meet the thermodynamic restriction (46a) is to model the flux of receptors by Fickian-diffusion, that linearly correlates h R R to the gradient of its chemical potential μ R R : by means of a positive definite mobility tensor M R . The following isotropic non linear specialization for the mobility tensor M R is chosen [60] where c max R R is the saturation limit for receptors, and u | R > 0 is the mobility of receptors. Definition (48) represents the physical requirement that both the pure (c R R = 0) and the saturated (c R R = c max R R ) phases have vanishing mobilities. Neither the mobility u | R nor the saturation concentration c max R R are assumed to change in time. Whereby experimental data indicate an influence of temperature, stresses, or concentrations, such a limitation can be removed without altering the conceptual picture. Noting that Fick's Law (47) specializes as where D | R = u | R R T is the receptor diffusivity.
The chemical kinetics of reaction (17) is modeled via the law of mass action (28). Experimental evidences [17] show that: (i) the equilibrium constant (22) is high, thus favoring the formation of ligand-receptor complex and the depletions of receptors and ligands; (ii) the diffusion of receptors on the cell membrane is much slower than interaction kinetics. Accordingly, it can be assumed that the reaction kinetics is infinitely fast, in the sense that the time required to reach chemical equilibrium is orders of magnitude smaller than the time-scale of other processes. For these reasons we assume that the concentrations of species are ruled by thermodynamic equilibrium at all times, and the concentration of complex c C R is related to the others by the equation (32a). This very same equation could be derived imposing Simple algebra allows deriving Eq. (32a), provided that to the equilibrium constant K (17) eq the alternative definition R is the standard Gibbs free energy.

Thermo-chemo-mechanics of cells
Cells show two main paradigmatic mechanical attitudes: active and passive. Active response is related to the ability of the cell to change, as a result of external cues, its own cytoskeletal conformation, i.e. to reorganize the morphology of the biopolymers net that provides the structural resistance during adhesion (to the ECM or to other cells), migration [61] (e.g. chemotaxis [62], mechanotaxis, and durotaxis [63]) and division (eg. mitosis). Passive, instead, refers to the mechanical response that each component of the cell has inasmuch material bodies, in accordance with their own internal structure and as a result of external actions.
Following again [29], the thermo-chemo-mechanics of endothelial cells can be stated stemming from energy and entropy balances. The referential specific Helmholtz free energy per unit volume ψ R T , c G R , c F R , C, ξ is taken as a function of temperature, strains (either C or E), concentrations c G R , c F R , and of some kinematic internal variables ξ that compare with the usual meaning in inelastic constitutive laws [6,22,[64][65][66][67]. Standard arguments-the so called Coleman Noll procedure-finally allow identifying the following thermodynamic restrictions The internal force, conjugate to ξ , will be denoted with the symbol χ , i.e.
Equation (51a) yields to the Clausius-Plank inequality, which under the assumptions of Curie symmetry principle [27] and thermal equilibrium, can be written as The flow of actin monomers is linearly related to the gradient of their chemical potential by Fick's assumption, consistently with the thermodynamic restriction (52b): The following positive definite, isotropic non linear specialization for the mobility tensor M G R is chosen [60] where c max G R is the saturation limit for receptors, and u | G R > 0 is the mobility of actin monomers. Assuming that the trapped species F has vanishing mobility is an alternative view of modeling the absence of their flux.

Decompositions.
The stress field S will be additively decomposed in the sum of the active and passive contributions, analogously to generalized Maxwell models Active response is related to cytoskeletal reorganization in stress fibers and pseudopodia, whereas the passive response reflects the mechanical behavior that each component of the cell has inasmuch material bodies. We base the theory for pseudopodia on a multiplicative decomposition of the deformation gradient Tensor F c , named swelling distortion is the local distortion of the material neighborhood of a point due to a volumetric swelling (de-swelling) due to the phase change of actin, from monomeric to a network of filaments and vice-versa. Its representation will be taken as F c = λ c 1, assuming therefore that a dense network of actin filaments form in pseudopodia. This approach conforms well for lamellipodia filament networks, although it might result inappropriate for slender and highly oriented microstructures seen in filopodia, which might be better captured by the protrusion-contraction uniaxial tensors presented in Refs. [54,68] or [69]. The following identities can be easily assessed: We assume that changes in J c occur because of changes in filaments J c = J c (c F R ) and define the partial molar volume of the pseudopodia as and it holdṡ The decomposition (56) leads to a multiplicative decomposition for the left Cauchy-Green tensor, too: with the swelling factor C c = J c2/3 1 and the elastic factor with a constant partial molar volume Ω C > 0.
In the realm of viscoelasticity, it is also common to perform a multiplicative decomposition of the deformation gradient F e into volumetric F e v and isochoric F e i factors The volumetric factor F e v = J e1/3 1 turns out to be completely identified by the determinant of F e , whereas the isochoric factor F e i = J e−1/3 F e obeys to the constraint det F e i = 1. The decomposition (62) leads to a multiplicative decomposition for the left Cauchy-Green tensor, too: with volumetric factor C e v = J e2/3 1 and the isochoric factor C e i = J e−2/3 C e .

Constitutive theory
The Helmholtz free energy density ψ R is modeled by decomposing it into separate parts: a thermal contribution ψ th R , a diffusive contribution ψ di f f R , an elastic contribution ψ el R , and an inelastic (also called configurational ) counterpart ψ in This splitting is here taken for granted without motivation. We will not indulge in the description of ψ th R (see Ref. [29] in case of interest) and we'll rather focus on the remaining parts.
Statistical mechanics depicts the entropy for isolated systems in terms of the density of states, the number of possible molecular configurations [70]. Making recourse to Stirling's approximation and since the entropy transforms with the volume by means of J , one finds that the following well-known expression of the entropy of mixing in the reference configuration arises: the universal gas constant R being the product of Boltzmann constant k B and Avogadro's number and having denoted with β = G, F and with ϑ β R the ratio We argued in Eq. (18) that, in view of the structure of the lipid membranes, the maximum number of moles or molecules per unit area for any species remains unchanged in time in the current configuration. The same argument does not seem to apply for the bulk, hence we take henceforth that is constant and write the free energy density for the continuum approximation of mixing [70] as (68) Note that if the saturation is constant in the current configuration, an explicit coupling of the free energy of mixing with the deformation arises by means of J . A new stress would come out, in view of the thermodynamic prescription (51a).
Following [64], we will define visco-elastic materials based on the multiplicative decomposition (63). Specifically, the free energy for visco-elastic materials will be defined as follows with ψ in R depending upon E e by means of C e i . The volumetric part of the elastic free energy is defined through J e , highlighting the role of the swelling tensor and of the concentration of pseudopodia, since in view of Eq. (61). On the other end, it holds hence C e i depends merely upon the state of deformation and not upon the concentration of species. This outcome reverberates upon the energetic contributions ψ el,iso R and ψ in R . The latter is such that a property physically grounded in the rheological model of Maxwell, for which we refer to [64] or [71]. Provided that the above holds, the selection for ψ el R and ψ in R is arbitrary. Their selection shall be different in modeling the passive behavior or the active response of pseudopodia and stress fibers. The elastic, reversible behavior that occurs once the viscous effects vanish (ideally at t → ∞ ) is captured by ψ el R . The inelastic free energy accounts for the non-equilibrium response due to viscosity-the so called dissipation potential. By thermodynamic restrictions (51) and identity (72) According to Eq. (73b), tensorial internal forces χ R can be interpreted as a non-equilibrium stress tensor of second Piola-Kirchoff kind, that accounts for the viscous response.
Inelastic internal entropy production (52a) was described by the internal flux variables ξ and by their energy-conjugate forces χ R . A simple way to satisfy constraint (52a) is choosing a positive definite operator L such that In case of isotropy, the fourth order operator L restricts to the scalar viscosity ν times the identity operator. Equations (73a), (74) provide evolution equations for χ R that allow the algorithmic integration of the constitutive law once a selection for the free energy densities ψ el R and ψ in R is made. The chemical potential of G-actin monomers and of Factin networks descends from thermodynamic prescriptions (51a), in the form While the chemical potential of actin monomers has merely an entropic nature, mechanical contributions enter the definition of the chemical potential of actin networks. Specifically, mechanics affects μ F R in the volumetric contribution ψ el,vol R through the swelling tensor C e v (70), whereas the isochoric tensor C e i was proven to be independent upon the concentration of species in Eq. (71). Nonetheless, the parameters of the viscoelastic loading-unloading law are expected to depend upon the extent of the polymerization reaction by means of the network concentration c F R in all terms of the mechanical free energy. The mechanical effect on the chemical potential does not propagate into the mass flux because the mobility of actin network is assumed to be negligible. Mechanics however enters the affinity of polymerization reaction (34) : the stress state is expected to favor polymerization nearby the lipid membrane and depolymerization towards the nucleus.

Multiscale scenario of cell viscoelasticity
Although the free energy scenario is rather clear, a specialization of the constitutive equations has not been attempted here and in many cases (as for the lamellipodia filament network and microtubules) it has not been attempted in the literature, to the best of our knowledge. The hindrance stands in the multiscale scenario of cell viscoelasticity: while the mechanical behavior and properties of single intermediate filaments, actin filaments, and microtubules has been nowadays quite clarified, at least in terms of relative stiffness and Fig. 4 Schematic picture of some cytoskeletal biopolymers strengths, bundles of the filaments, their response, polymerization, shape and time evolution is hard to be captured by comprehensive models at the "macroscopic" scale through appropriate free energies. As a consequence, the ability of models to capture the mechanics of fundamental cellular processes (as chemotaxis, cell sprouting, junction and differentiation [72], endocytosis [73] and exocytosis to cite a few) still requires abundant research before gaining predicting capabilities in simulations.
The cytoskeleton, an interconnected network of regulatory proteins and filamentous biological polymers depicted schematically in Fig. 4, undergoes massive reorganization during cell deformation, especially after cell rolling and adhesion [37,74] and in mediating, sensing and transduction of mechanical cues from the micro-environment [75]. Homogenized models for the mechanical response of a cell shall condense in effective properties the: i) polymerisation/depolimerisation of filaments; ii) the process of cross-linking that determine the architecture of cytoskeletal filaments; iii) the passive mechanical properties of the cytosol. The thermodynamics of statistically-based continuum theories for polymers with transient networks [76][77][78][79][80] appear to be a valuable candidate for the selection of free energies ψ el R (c F R , C) and ψ in R (c F R , E, ξ ). At present however, such a comprehensive model has not yet been proposed for the pseudopodia driven cell motion. Classical models as hyperelastic Saint-Venant [54] or newtonian viscous fluids [81] eventually surrounded by a hyperelastic, zero-thickness membrane [82] have been used for the pseudopodia, whereas a very large amount of literature concerns pseudopod dynamics ( see for instance [83] and the large literature therein ) or ameboid motion [84]. Different approaches to cell motility, as for active gel theory coupled to the classical theory of thin elastic shells, are also widely used [85], but are not discussed in this work. The framework described herein, including myosin dynamics, phase transformations between G-actin and F-actin, has been depicted in a set of publications by the group of H. Gomez [36,86]. The flow of the F-actin network was treated as a Newtonian fluid and directed by its velocity. A one dimensional yet comprehensive model has been proposed in [87].
The multiscale scenario is invoked also for cell contractility. There are evidences [88] that the interaction among filaments, motors, and cross-linkers is mechanically stimulated. As reported in [75], myosin binding to actin fibers occurs in a force-dependent manner, as well as the contractile response of actomyosin to extracellular stiffness. According to [89], force feedback controls motor activity and increases density and mechanical efficiency of selfassembling branched actin networks, thus suggesting that those feedbacks could allow migratory cells adjusting their viscoelastic properties to favor migration. Mass transport and cell contractility have been accounted for in several publications with different degree of complexity [50,69,79]: to the best of our knowledge, however, the force transmission has always been modeled stemming from the similarity between the sarcomeric structure of stress fibers and the actin-myosin interactions in muscle cells. In [35] a multi-dimensional network of stress fibers was built on the notion of a representative volume element, in which stress fibers can form in any direction with equal probability. An average macroscopic stress is then recovered from the fiber tension, which in turn is generated by the cross-bridging cycles and described by a Hill-like relation [90] of viscoelastic nature. Experimental evidences, however, seem to show that such a resemblance might be questionable in the dynamics and mechanics of endothelial cell spreading [19] and hence that the predictive capability of this family of models might be poor for this family of cells.
Finally, the passive response of the cytosol, provided mainly by the intermediate filaments attached to the nuclear and plasma membranes, has been modeled by several authors by means of classical models as linear elasticity [79], the finite strain generalization of Hooke's law [35] or a Neo-Hookean potential energy where G 0 is the initial shear modulus and G ∞ is the shear modulus at the end of the viscous processes. This classical choice of Helmholtz free energy is associated to efficient integration schemes, depicted in [71].

Application
Although the focus of the present note stands in the establishment of a framework, rather than providing a model for any specific phenomenon, it seems of use introducing one simple but concrete example of a specific problem, in which the novelties introduced in the previous section can be used. Specifically, we will briefly describe the relocation of vascular endothelial growth factor receptors-2 (VEGFR2) observed in an in vitro experimental setup and replicated numerically. Such a problem was discussed in [17,28] in the framework of small displacements and strains and will be extensively presented in a companion paper [91] using the finite strain multiphysics framework developed herein. VEGFR2 is a pro-angiogenic receptor expressed on endothelial cells (ECs) and is the main mediator of the angiogenic response. The interaction between VEGFR2 and extracellular ligands, produced by tumor cells, is essential to cancer growth. Specifically, ligand stimulation causes the relocation of VEGFR2 in the basal aspect in cells plated on ligand enriched extracellular matrix both in vitro and in vivo, and ultimately receptors-ligands interactions activate the ECs division and proliferation towards tumor cells. Upon release, growth factors associate with the extracellular matrix and act as ECs guidance during neo-vessel formation.
The binding-unbinding interaction (17) was tailored to describe the interaction between VEGFR2 (R) and VEGF (L), which produces a receptor-ligand complex (C). Mass balance equations along the advecting membrane lead to the chemo-transport problem (29), which comprises the governing equations for the relocation on the membrane. The geometrical evolution of the cell (and hence of the membrane as well, onto which VEGFR2 relocation takes place) is governed by balance equations (44) coupled to the thermodynamic restrictions elucidated in Sect. 6. At present, research on the identification of suitable free energies ψ el,iso for the active stress S active is on going, based on statistical mechanics of cytoskeletal reorganization in stress fibers and pseudopodia. We thus accounted only for the passive stress S passive by means of rubber visco-elasticity (76). Such a multi-physics initial boundary value problem in the bulk and on the membrane of the cell, rephrased in a weak form and further discretized via finite elements, has been implemented in a high performance computing code with a staggered Newton-Raphson solver, in the deal.ii framework ( http://www.dealii.org ).
The resulting code has been used to simulate the relocation of VEGFR2 expressed on endothelial cells on the cell membrane during the mechanical adhesion and spreading of cells onto a ligand-enriched substrate. A co-designed experimental and computational campaign unveiled the multiphysics progression of the process. The geometrical evolution of the cell was recorded for 2 hours in time-lapse microscopy. During this timespan, three mechanically relevant events could be identified: the floating and adhesion of the EC on the ligand rich μslide, and eventually the spreading onto the latter. Three limiting processes characterize the depletion of VEGFR2. In the first phase, Fig. 5(b)-(d), the VEGFR2 depletion is dominated by the chemical interaction between receptor and ligands during the adhesion between the cell membrane and the ligand-reach substrate. The second phase, up to 600 s, illustrated in Fig. 5(d)-(i), is driven by the mechanical spreading of the cell: the cell-substrate contact dynamics stimulates the formations of complexes, since the mechanical spreading makes new free receptors available for the binding with the ligands. Eventually, the last phase is dominated by a lower complex formation rate (from 600 s to 7200 s) and takes place after the EC spreading was completed, thus resulting transport-dominated, see Fig. 5(i)-(j). Free VEGFR2, guided by concentration gradients, move from the apical part of the cell towards the basal one, where the binding-unbinding interaction (17) occurs. Ultimately receptors depletion is complete on the cell membrane-see Fig. 5(j) : such an event is actually unrealistic, since an immobile fraction of VEGFR2 shall be accounted for, as it will be illustrated further in [91].

Concluding remarks
In this note, a multi-physics framework of protein relocation on the advecting lipid membrane during cells spreading and motion has been put forward. It sets the (continuum) thermodynamic background for simulations of receptor recruitment during migration: simulations carried out in [17] stem from a simplified form of the framework and described the limiting factors in vascular endothelial growth factor receptors relocation; similarly, we discussed in [18] the relocation of integrins on the membrane and their interactions with growth factor receptors; a companion paper [91] deals with the relocation of vascular endothelial growth factor receptors on advecting lipid membrane during endothelial cell adhesion and spreading. Those simulations may have a significant impact in biology and in the pharmacological treatment of cancer, either in view of their predictive nature in virtual experiments, or by clearly identifying the sequence of processes that limit the relocation of targeted proteins during in vitro experiments.
The present work still has significant limitations, yet by illustrating a complex and rigorous scenario it might be a cornerstone to account for several further processes. To cite a major phenomenon that has been insufficiently discussed here, the proteins transport on the membrane is crucially coupled to the cytoskeleton reorganization, which is related to the motion of integrins on the membrane: the formation of focal adhesion sites is preliminary to stress fibers generation and contractility. Internalization of complexes is another occurence not included in this work. Further publications, therefore, will be devoted to extend this framework to these and others challenging tasks.
We also aimed in this paper at recollecting recent publications from several schools on cell mechanics, encasing them in a unified framework, being aware that a comprehensive account of publications is significantly hard in view of the broadness of the literature in the field. We clarified that for some processes, as for contractility and protrusion, either a thermodynamically consistent formulation has not been devised yet or it stems upon simplistic models that do not account for the microstructural evolution of the biopolymers. Even in this fascinating field, the last word is far from being spoken.