An elasto-plastic biphasic model of the compression of multicellular aggregates: the influence of fluid on stress and deformation

We present a mathematical model of the compression of multicellular aggregates, and we specialise it to a compression-release test that is well known in the biological literature. Within the adopted mechanical setting, a multicellular aggregate is studied as a biphasic system consisting of a soft solid porous medium saturated with an interstitial fluid. In particular, together with the deformation of the considered aggregate, the characterisation of the model outlined in this work relies on four fundamental features. First, by assuming the interstitial fluid to be macroscopically inviscid and to evolve according to the Darcian regime, we resolve its flow and determine the associated time dependent pressure distribution. Second, we focus our attention on the remodelling of the compressed aggregate, that is, on the rearrangement of its internal structure in response to the external loads applied to it. Specifically, we look at the way in which such a rearrangement is induced by the considered experiment and at how it affects the mechanical behaviour of the aggregate. Moreover, we introduce a remodelling-dependent permeability tensor with the purpose of visualising a more direct influence of remodelling on the dynamics of the aggregate’s interstitial fluid. Finally, we resolve the interactions exchanged between the aggregate and the compressive apparatus. This task necessitates the formulation of an appropriate contact problem, thereby calling for the description of the evolution of the area through which the aggregate and the apparatus exchange mechanical interactions. In particular, the continuity conditions to be applied on such a contact area are discussed. Our numerical simulations show the role played by the different phenomena accounted for in the model and the overall dynamics of the aggregate within the considered experiment.


Introduction
The way in which living systems evolve in time, interact with their surrounding environment and sustain their biological functions is the result of processes of different nature, some examples of which are given in [35,36,38,72,78]. Because of the difficulty in keeping track of all these processes, to characterise biological structures such as tissues or organs, it is important to conceive benchmark experiments on systems that, although comparatively simpler, permit to observe the most relevant aspects of a given phenomenology. Within the context of the biomechanical research on cancer, in which a prominent role is given to the mechanical aspects of tumours in the early stages of their formation, multicellular aggregates are largely employed for the purposes of benchmarking described above [1,22,25,56,[61][62][63][97][98][99]. Indeed, they can be obtained by culturing cells in vitro [22,86] and are commonly referred to as multicellular spheroids in the literature (see e.g. [1,13,22,25,26,39,69,72,78,86,99], since they are given a spheroidal shape, when they are prepared for experiments. Their relatively simple geometry and internal structure allow for mathematical models that capture their most fundamental mechanical features, but that are free of the mathematical complications necessary for other systems [57,61,68,69,73]. More specifically, multicellular aggregates/spheroids are a good compromise between one-dimensional models of tissues, as is the case of Fig. 1. A representation of the multicellular aggregate (before the compression begins) and of the plates. The sphere represents the multicellular aggregate prior to the commencement of the compression, and the cylinders represent the two plates of the experimental apparatus. The aggregate finds itself inside the inner chamber of the experimental apparatus, whose walls have not been drawn in the picture (see e.g. [34]) at a certain fixed temperature by circulating water, pumped from a system located in the outer chamber [32][33][34]55,92]. A representation of the aggregate and of the plates is shown in Fig. 1 (for the visualisation of the experimental apparatus, the Reader is referred to, for example, [34]).
We remark that, in the considered experimental framework, the multicellular aggregate is prepared in such a way that its shape is spheroidal prior to the onset of deformation and, in fact, in the simulations presented below, it is described as a sphere before the deformation commences. This is why the aggregate is commonly referred to as "spheroid" in the biological literature. However, even before starting the compression, spherical symmetry is broken by the force of gravity.
The compression starts after the multicellular aggregate is placed inside the inner chamber of the experimental apparatus and put in contact with the upper and with the lower plate [32][33][34]55,92]. At this initial stage, only the north and south pole of the aggregate are put in contact with the plates and, specifically, in correspondence to the centre of the lower face of the upper plate and to the upper face of the lower plate, respectively [32][33][34]55,92]. Then, the experiment is conducted in displacement control: the upper plate is kept fixed, and a displacement ramp, defined as a function of time, is applied to the lower plate. This results in a compression of the aggregate between the two plates, which is estimated in terms of the force acting between the compressed aggregate and the upper plate [32][33][34]55,92]. Depending on the applied displacement ramp, different compressive profiles of the multicellular aggregate are possible, and this is useful to understand more deeply how it responds to external loads and to better estimate its mechanical properties. In addition, the procedure permits to observe the possible recovery of the aggregate's initial shape when the imposed ramp is removed [32][33][34]55,92].

Modelling hypotheses
As motivated in [42], we adhere to a Continuum Mechanics approach to investigate the mechanical response of the considered multicellular aggregate, and we summarise the main reasons for which describing aggregates as continua is appropriate for the class of problems addressed in this work. In fact, the aggregates we refer to in our model are characterised by a radius that, approximately, ranges from "hundreds of micrometres up to some millimetres" [42], as is the case of tumour spheroids in avascular stage [32][33][34]55,92], whereas the cell radius is about 5 μm up to approximatively 10 μm. This implies that, in general, an aggregate of such a size comprises thousands or hundreds of thousands of cells [32][33][34]42,55,92]. Consequently, the employment of discrete models may become too demanding. Moreover, the radius of the aggregate, as well as the magnitude of the contraction to which it is subjected, is characterised by length scales that are sharply separated from the typical length scales at which intracellular, cell-cell and cell-ECM interactions occur. In light of all these considerations, and for the physical ZAMP An elasto-plastic biphasic model of the compression Page 5 of 39 79 phenomena and experimental procedures that we are taking into account, a continuum approach seems justified [32][33][34]42,55,92].
As is the case in many studies regarding, for instance, the mechanical characterisation of tumours or soft biological tissues, we resort to the theory of binary mixtures [6,10,15,41,45], according to which an aggregate is assumed to comprise a solid and a fluid phase. The solid phase is mainly constituted by the cytoskeleton and other components, such as the ECM, the nuclei of the cells and the other individual sub-cellular units, which find themselves in the cytosol. The fluid phase is constituted by the interstitial fluid present in the pores of the ECM and among the cells [6,10,15,41,45] as well as, in principle, by the fluid in the cells. However, as anticipated in Introduction, we emphasise that there are at least two main reasons for giving no prominent role to the intracellular fluid in our model. The first reason is that we do not regard its dynamics as distinguishable from the dynamics of the cells as a whole. This is related to the fact that, in our approach, the intracellular fluid is one of the constituents of what we call "solid phase" (in fact this should be called "effective solid phase") and contributes to determine its overall mechanical properties. In passing, we notice that this formulation is coherent with the one developed in [12] for multi-constituent, biphasic mixtures that are used to describe fluid-saturated porous media through the so-called Hybrid Mixture Theory. More in detail, besides the intracellular fluid, in our model the solid phase comprises also collagen as well as the structural parts of the cells (e.g. membrane, nucleus, organelles, endoplasmatic reticulum, etc.) and of the ECM, and these constituents are all constrained to have the same velocity. The latter coincides with the velocity of the centre of mass of the solid phase of the aggregate and is kinematically distinguishable from the velocity of the extracellular fluid, whose dynamics is studied in the sequel under the assumption of validity of Darcy's law. The second reason is that we study no mass exchange between the extracellular fluid and the intracellular one. In fact, such exchange of mass could contribute to the mechanical response of the aggregate in the proximity of the contact zones with the compressive apparatus. However, we regard this phenomenon as secondary since, for the time being, we focus our attention on what can be predicted by considering the interstitial fluid only. On the contrary, processes of this kind must be taken into account, for instance, in chemo-mechanical problems that involve the transfer of mass between the fluid and the solid phase of a cellular aggregate, or when studying growth [5,69]. In the latter case, however, the time scale of growth is of several days (e.g. from 1 day to 20 days), whereas the time scale of the experiments addressed in this work is of the order of minutes.
To describe mathematically the type of systems depicted above, we make the following hypotheses. First, we assume that the mechanical behaviour of the collagen, of the ECM and of the structural parts of the cells results macroscopically into a hyperelastic behaviour of the solid phase of the aggregate. Second, we hypothesise that the interstitial fluid is macroscopically inviscid and, as anticipated above, that it obeys Darcy's law. Moreover, in spite of its importance in the overall description of aggregates, in this work we neglect the dynamics of the macromolecules transported by the fluid so that some chemical processes are not resolved. This, in turn, restricts our present model to those situations in which chemical phenomena are weakly coupled with the mechanical ones. Note, however, that this limitation is acceptable for the compression-release experiment that we are describing (see e.g. [32][33][34]), which, indeed, looks only at merely mechanical aspects of the aggregate related to the imposed loading history.
The main focus of this work is the remodelling of the considered multicellular aggregate in the context of the simulated experimental procedure. We recall that, with the term remodelling, we refer to a family of transformations that take place in a biological system and lead to the adaptation of its internal structure to both external and internal stimuli of different nature [20]. Hence, in conjunction with the deformation of the aggregate, we also look at how the simulated experiment yields structural rearrangements of the aggregate at the cell scale [32,34,91,93,100]. Note, however, that no phase transitions are considered for the time being.
Remodelling is typically associated with anelastic distortions, as is also the case for other anelastic processes, such as volumetric growth or mass exchange among the constituents of a tissue [4,41,45,48,69].
In the case of a multicellular aggregate, remodelling manifests itself through inter-cellular and intracellular interactions [32,40,48]. The latter ones encompass, for instance, the transformation of the actin filament network within the cells, the disaggregation and formation of adhesion bonds among the cells and between the cells and the ECM, as well as the change of shape, position and orientation of the cells [40,48,97]. All these phenomena, occurring at the cellular scale (or even at finer length scales), produce the mechanical response of the aggregate to the applied loads, in terms of stress and of the structural evolution associated with it.
We assume that the type of remodelling occurring in the performed experiment does not produce any variation of mass [2,32,55,66,82,96]. In addition, its characteristic time scale is set by the time scale of the applied loads and, more generally, of the experimental procedure [2,32,55,66,82,96]. For the study at hand, and by considering the experimental protocol outlined in Sect. 2.1, the remodelling of the multicellular aggregate can be conceived to occur over a time scale of the order of seconds or minutes [2,32,55,66,82,96]. Hence, such kind of structural reorganisation is well separated from the growth of the aggregate itself, which, in general, takes place over days or weeks [2,32,55,66,82,96]. In fact, this is the time scale at which the mitosis and the apoptosis of cells occur [2,32,55,66,82,96]. Hence, we can study remodelling by neglecting growth, both volumetric and appositional, and by decoupling the anelastic distortions associated with remodelling from those pertaining to growth [2,32,55,66,82,96].
As for growth, remodelling can often be decoupled, with analogous motivations, also for other biological processes, such as the mass exchange between the solid and the fluid phase, cell differentiation and morphogenesis [2,32,55,66,82,96]. We remark, however, that this is a modelling hypothesis, strictly related to the specific phenomena and experiments studied in this work.
Starting from the evidence that the different sub-cellular units have properties that make them substantially dissimilar one from the others, the mechanical properties of a single cell are distributed inhomogeneously within it [40,98]. Consequently, also the overall mechanical behaviour of an aggregate is, in principle, inhomogeneous. However, since we do not know, up to now, how the elastic parameters depend on the material points and, possibly, on time, we assume them to be constant. This modelling choice has been tested in previous works [40,42], where the numerical results are in good agreement with the experimental ones. Moreover, for the kind of aggregates that we are studying, the amount of the ECM is low with respect to the cellular component and to the quantity of cross-links among the cells [43]. This experimental evidence shows that the possibility for collagen reinforcing fibres to anchor themselves to the ECM is limited [43]. In this respect, we can assume that the multicellular aggregates addressed in this work display isotropic mechanical response.
Finally, as done in [42], we assume the compressive plates to be linear elastic, isotropic, homogeneous, consisting of the same non-adhering material [32][33][34]55,92], and sharing a slightly permeable, time dependent contact area with the multicellular aggregate. In addition, we require the boundaries of the compressive plates that are not in contact with the aggregate to be impermeable. We emphasise that, in this work, we model the plates as deformable saturated porous media, rather than as rigid solid bodies. This choice has been dictated partially by technical reasons, as will be explained in Sect. 7.3, and partially by the need of a formulation that can be easily employed to a variety of experimental set-ups.

Topology of the system
In the sequel, with the wording mechanical system we refer to the complex consisting of the multicellular aggregate and the two plates of the compressive apparatus, and we observe its evolution in the threedimensional Euclidean space S , over the interval of time T = [0, T ], for some T > 0.  2. A cross section passing though the symmetry axis of the mechanical system at time t = 0 s and at a generic time t ∈ ]0, tramp[ of the compression ramp. Note the relations Σ (0) = I (0) and Σu(0) = Iu(0) to indicate that at t = 0 s, there is no mechanical contact (in the sense explained later, i.e., absence of mechanical interactions) between the aggregate and the upper plate, while, because of gravity, there is mechanical contact between the aggregate and the lower plate. The red dots indicate the view of ∂∂ T Bu(0) and ∂∂ T B (0) in the cross section. In the picture on the right, the boundaries ∂ T Bu(t), have been omitted for the sake of a lighter representation

Current configuration
We denote by B(t) ⊂ S the portion of S occupied by the upper plate, by the multicellular aggregate and by the lower plate at time t, i.e., where B u (t), B a (t) and B (t) are open and connected subsets of S . Note that, in principle, these three regions may be detached from one another over a subset of the time window T . However, for the case studied in this work, the boundaries of B (t) and B a (t), denoted by ∂B (t) and ∂B a (t), respectively, have at least one point in common for any t ∈ T . Hence, it holds that and, in the course of the experiment, I (t) may either degenerate in a point or evolve into a surface. On the other hand, for the way in which the experiment is prepared (the aggregate and the upper plate touch each other prior to compression), if we indicate with ∂B u (t) the boundary of B u (t), the intersection is empty in part of the release phase and in the relaxation phase of the load, while it is non-empty otherwise, and it may either evolve into a surface or reduce to a point. Because of these facts, we define the configuration of the mechanical system at time t ∈ T as A graphical description of the initial and of the current configuration of the system is given in Fig. 2. We emphasise that Eqs. (2) and (3) rely exclusively on the topological characterisation of C (t), which is related to the evidence according to which B a (t) shares some portions of its boundary with those of B (t) and B u (t), and we may refer to these situations as topological contact. Yet, we do not speak of mechanical contact unless mechanical actions occur through I (t) and I u (t). Conversely, when I (t) and I u (t) are non-empty and host mechanical interactions among the aggregate and the plates, we do speak of mechanical contact zones or mechanical contact areas (from here on, just contact zones or contact areas), and say that B a (t) is in contact with B (t) and/or with B u (t). Moreover, to emphasise that contact requires both a non-empty intersection between the boundary of the aggregate and that of at least one of the plates as well as the exchange of mechanical actions through such intersection, we denote the contact zones by Σ u (t) and Σ (t). Hence, if I (t) and I u (t) are non-empty and feature mechanical interactions, we set However, when I β (t), with β ∈ {u, }, is empty, or no mechanical actions take place on it, we conventionally set Σ β (t) equal to the empty set. We remark that in our simulated experiment, Σ (t) is always non-empty, because I (t) is always non-empty and hosts mechanical interactions at each instant of time t ∈ T . Differently, Σ u (t) is empty in part of the release phase and in the relaxation phase of the load, i.e., for those instants of time for which I u (t) is empty, and prior to the application of the external load, i.e., when I u (t) consists of a point but no mechanical interactions occur on it. The contact areas in Eqs. (5) change in time accordingly to the loading history, deformation, fluid flow, distortions associated with remodelling and material properties of B a (t), B u (t) and B (t). In particular, Σ (t) and, when non-empty, Σ u (t) could be a singleton or a surface. When Σ (t) and Σ u (t) are singletons, the contact of B a (t) with B (t) and B u (t) occurs in correspondence of the south pole of the aggregate and the centre of the upper face of B (t) and of the north pole of the aggregate and the centre of the lower face of B u (t), respectively. Finally, since Σ u (t) and Σ (t) are disjoint subsets of ∂B a (t) for all t ∈ T , the overall contact set can be defined as which, depending on the instant of time t ∈ T , could be one surface only, the disjoint union of two surfaces, the disjoint union of a surface (i.e., Σ (t)) and a singleton (i.e., Σ u (t)), or a singleton only. In the sequel, we exclude the case in which Σ(t) is the union of two singletons because, in such a case, even though I (t) and I u (t) are simultaneously singletons, no mechanical interactions occur, since this situation is prior to commencement of the experiment. We remark that the topology of C (t) is strongly influenced by the circumstance for which I u (t) may be the empty set. Indeed, if I u (t) = ∅ (so that Σ u (t) is empty too, and Σ(t) = Σ (t)), C (t) consists of the disjoint union of two connected sets, which are B u (t) and B a (t) ∪ Σ (t) ∪ B (t). On the other hand, for I u (t) = ∅, C (t) is a single connected set, comprising the two plates, the aggregate and the contact set Σ(t), i.e., For further use, we also introduce the portions of the boundary as the portions of ∂B u (t) and ∂B (t) that are instantaneously not in contact with the aggregate. Analogously, we set Γ a (t) := ∂B a (t)\Σ(t) (8) to indicate the portion of ∂B a (t) deprived of the contact area.

Reference and initial configuration
As done in the literature from which we took inspiration for this work [40,42,46,79,80], we claim that a reference configuration exists for the mechanical system at hand. As well-established in Continuum Mechanics, this hypothesis is usually made because one can solve the equations of a given problem on a known, fixed geometry. In our case, there is also the advantage that one can describe the evolution of the aggregate's internal structure by having recourse to the standard theory of elasto-plasticity. From the computational point of view and, specifically, according to the requirements of the software used for our simulations, the assignment of the reference configuration for the aggregate and the plates amounts to supplying the reference geometry from which the initial configuration is determined. The latter, in turn, is identified with the configuration of the system prior to the commencement of the imposed displacement of the lower plate and is determined by solving the model equations. Note that, this is, in fact, a sort of initialisation of the simulation of the true experiment. Finally, for the problem considered in this work, the reference configuration is assumed to be undeformed and stress-free, while the initial configuration features a non-trivial stress distribution and a small amount of elastic deformations. This is due to gravity, which, in contrast to what is usually done for relatively soft tissues, is not neglected in our model. On the other hand, the stress associated with the initial configuration is not sufficient for plastic-like distortions to appear within the aggregate, so that no evolution of its internal structure is set off. For the reasons outlined above, although the multicellular aggregate has the shape of a sphere of radius R in its reference configuration (see also Table 1), it loses the spherical symmetry already in the initial configuration. Consequently, we prefer to speak of "multicellular aggregate" rather than of "spheroid" in the description of its evolution. Clearly, the deviation from the spherical shape due to gravity depends on the material parameters of the multicellular aggregate and may be negligible in some cases.
It is worth to mention that the introduction of a reference configuration for a biological system undergoing remodelling is supported by the literature (see, for instance, [24,40,42,46,79,80]). This is particularly true when a biological medium is not subjected to appositional growth, as is the case for the multicellular aggregate studied in this work.
Similarly to what we have done in Sect. 3.1, we denote by B ⊂ S the union where B u , B a , and B are the reference configurations of the upper plate, aggregate and lower plate, respectively. Moreover, given ∂B u , ∂B a and ∂B , the intersections are non-empty, and each of them corresponds to a singleton. We remark that, at this stage, no mechanical actions are supposed to be exchanged through I or I u , so that we cannot speak of "contact areas". Finally, we note that the set C := B ∪ I u ∪ I is connected and we naturally identify it with the reference configuration of the system as a whole. Given t 0 = 0 s as the instant of time at which the experiment begins, i.e., the time at which the imposed displacement starts, we introduce the embedding χ 0 : B → S , so that one can identify B α (t 0 ) = χ 0 (B α ) with the initial configuration of B α , with α ∈ {u, a, }, and one can define Now, once the boundaries ∂B u (t 0 ), ∂B a (t 0 ) and ∂B (t 0 ) are introduced, the intersections are non-empty and correspond to a singleton and to a "small" contact zone, respectively. At this stage, indeed, the contact occurs only between the aggregate and the lower plate and, as reported above, such an exchange of mechanical actions is induced by the presence of the gravity in the model. Finally, we define the initial configuration of the overall system as the union C (t 0 ) := B(t 0 ) ∪ I u (t 0 ) ∪ I u (t 0 ).

Model equations
Given a point x ∈ S and a point X ∈ B, we denote by T x S and T * x S and by T X B and T * X B the tangent and the co-tangent spaces of S at x and of B at X, respectively [67]. Moreover, for all x ∈ S and for all X ∈ B, we introduce the metric tensors g(x ) :

Balance laws
We model the multicellular aggregate and the two plates as biphasic mixtures, comprising a porous solid matrix and an interstitial fluid. For a system of this kind, the local form of the mass balance laws of its phases reads where the superscripts "(s)" and "(f)" refer to the solid phase and to the fluid phase. In Eqs. with and are constrained by the saturation condition ϕ From here on, we assume that the intrinsic mass densities are constant (so that both phases are incompressible) and, by manipulating Eqs. (13) according to a rather standard procedure (see e.g. [5,51]), we can recast them as (no sum over α) where D α is the substantial derivative of ϕ (s) α with respect to the solid phase motion and q α := ϕ α ] is said to be the filtration velocity. Note that, although the aggregate undergoes remodelling, Eqs. (16a) and (16b) are equally valid for the aggregate itself and the two plates. This is because the type of remodelling considered in this work does not produce any variation of mass.
If no body forces other than the force of gravity act on the system, and if all the inertial terms of the system are negligible (including those associated with the acceleration of the lower plate), the balance of linear momentum for the solid and the fluid phase reads where σ (s) α and σ  By accounting for the fact that the mixtures in the aggregate and in the plates are closed with respect to linear momentum, so that the condition m (s) α +m (f) α = 0 applies, Eqs. (17a) and (17b) can be combined to obtain In the sequel, we adhere to a rather "classical" description of the mechanics of the biphasic mixtures that are commonly encountered in the study of biological problems [41,48,79]. For this reason, we approach the balance laws (16a), (16b), (18a) and (18b) under the following hypotheses, which apply both in the aggregate and in the plates: • The solid phase is hyperelastic. Hence, its Cauchy stress is σ is referred to as the constitutive part of σ (s) α and is obtained from a suitable strain energy density, as shown in the next section.
• The fluid evolves in the so-called Darcian regime. Accordingly, σ (f) α is hydrostatic and reads σ [45,46,48,51], where m with k α being an invertible second-order tensor referred to as permeability tensor. We remark that, in the present context, k α is assumed to be symmetric and positive definite. Note also that the expressions of σ These hypotheses allow to replace the balance laws (16a), (16b), (18a) and (18b) with We point out that, since Eq. (19a) expresses the evolution of ϕ (s) α exclusively in terms of the volume changes driven by v (s) α , it is, in fact, decoupled from the remaining two equations, which are, thus, those that have to be implemented. However, since remodelling occurs in the aggregate, this has to be done in conjunction with a suitable remodelling law, which will be introduced in the next sections.

Kinematics of large deformations and "material form" of the balance laws
For the considered system, we introduce the motion of the solid phase as the one-parameter family of embeddings [67] so that B(t) = χ(B, t), for all t ∈ T , and χ( · , t 0 ) = χ 0 . Moreover, upon defining the restrictions of (2) and (3) can be determined from χ as For future use, we also define the maps χ : B × T → S and χ α : B α × T → S , and we denote by the deformation gradient tensor associated with χ α , along with its determinant J α (X, t) = det F α (X, t).
We recall that J α satisfies the differential relationJ α ( so that the Piola transformation of the filtration velocity reads where Φ is the Piola transformation of the volumetric fraction of the fluid phase. Note also that the Piola transformation of the volumetric fraction of the solid phase is defined analogously, i.e., as Φ Finally, we also introduce the right Cauchy-Green deformation tensor, defined as the pull-back of the metric tensor g [67], i.e., with g α (X, t) = g(χ α (X, t)). By virtue of the definitions given so far, Eq. (19a) reduces to the ordinary differential equationΦ (s) α = 0, to be solved for all α ∈ {u, a, } and whose solution is Φ αR (X) is referred to as referential volumetric fraction of the solid phase and is, by definition, constant in time. Furthermore, the remaining balance laws (19b) and (19c) become, in material form, where p α is the pressure expressed as a function of material points and time, i.e., p α (X, t) = p α (χ α (X, t), t), and we denote the Piola transformation of the permeability tensor (see e.g. [11,21,31] and references therein) and the constitutive part of the solid phase first Piola-Kirchhoff stress tensor as Both quantities are supplied constitutively, as will be shown in the following sections.

Bilby-Kröner-Lee decomposition and law of remodelling
To provide a mathematical characterisation of the rearrangement of the aggregate's internal structure, we refer to the theory of elasto-plasticity and, by following the literature (see [9,59,65,71,89]), we invoke the Bilby-Kröner-Lee decomposition of the aggregate's deformation gradient tensor F a . In particular, we claim the existence of two non-integrable second-order tensors, denoted by F e and F p , such that where F e is the tensor of elastic distortions, also called "accomodating tensor", while F p is the remodelling tensor, which describes "plastic-like distortions" [21,24,48]. Both tensors are non-singular, and their determinants are strictly positive and denoted by J e and J p , with J p = 1 under the hypothesis of isochoric remodelling. Note that, in this case, the identity J a = J e J p = J e holds true. We emphasise that remodelling is the phenomenon, while the terminology plastic-like distortions indicates the description of the phenomenon itself, whose manifestation is here conceived through distortions that resemble, to a certain extent, the plastic distortions in elasto-plastic materials.
We remark that, for the sake of a lighter notation, we do not add the subscript "a" to F e and F p . Indeed, these two tensors are defined only in B a , and not in B u or in B , since remodelling characterises the aggregate only.
For completeness and future use, following the same notation as in [24], we introduce the vector space We recall that the elements of N X (t) are the elements of the tangent space T X B a as they look like when they relax at time t, and are therefore said to be in their natural, or ground, state [44]. Many further details can be given about the BKL decomposition [59,87], its use in Biomechanics [54,79,84], and the most recent advances in the comprehension of its connection with residual stresses [3,17,32,34,37,58,60,70,77,93,94].
Here, however, we do not discuss these topics, since they are not the focus of our present work.
One can construct the Cauchy-Green remodelling tensor with x = χ a (X, t) and X ∈ B a , so that b e is understood as the push-forward of B p through the motion χ a . At this stage, we can compute the Lie derivative of b e along the velocity v (s) a , i.e., Note that, for the sake of a lighter notation, from here on we simply set Lb e (x, t) ≡ L v (s) a b e (x, t), with x = χ a (X, t) and X ∈ B a . We recall that the Lie derivative of b e provides an objective rate of the variation of b e , evaluated along the integral curves of v (s) α [67]. In the sequel, the solid phase is assumed to be isotropic with respect to both its elastic and hydraulic properties. Consequently, the stress tensor σ (sc) α and the permeability k α can be expressed constitutively as functions of b e only [64]. Moreover, the distortions associated with the remodelling can be described by B p , when studied from the point of view of the reference configuration B a , and by b e , when studied in the current configuration B a (t). Indeed, B p measures the structural changes of the tissue in terms of the change of the material metric from G to C p , and, thus, from We remark that the results discussed above are peculiar of the constitutive context of our work because of the hypothesis of isotropy, but, for example, they are unsuitable in the case of anisotropic media or when remodelling occurs through inelastic rotations. Indeed, if the remodelling distortions turned out to be rotations, B p would be equal to G −1 and no change of metric would be observed. In fact, by constructing our model exclusively on B p , we lose information on plastic-like rotations, because they are unresolved in our study, as shown in the next section.
A last remark concerns Eq. (31), which, consistently with what has been said above, prescribes that the Lie derivative of b e is driven by the time change of B p , thereby indicating that the model requires an evolution law for this tensorial variable. Such an evolution law can be equivalently prescribed in one of the two formsḂ where Y and y are suitable generalised forces power-conjugated to B p and b e , respectively, and the ellipsis points stand for any collection of additional variables that modulate the evolution, such as the model parameters.

Constitutive framework
Because of the hypothesis of hyperelastic solid phase both in the plates and in the aggregate, the stress tensor P (sc) α can be derived from a Helmholtz free energy density, which has to be specified for each α ∈ {u, a, }. Moreover, we assume both the plates and the aggregate to be homogeneous, so that the material parameters involved in the definition of the Helmholtz free energy density can be taken as constants.
To prevent the occurrence of compaction, we endow both the plates and the aggregate with a strain energy penalty term. Since remodelling is assumed to be isochoric in this work, and, thus, it holds that J e = J a in the aggregate, such penalty term can be expressed as a function of J α , for α ∈ {u, a, }, and it can be given the same functional form for all such compartments. In fact, since preventing compaction amounts to enforcing that J α (X, t) may vary only in the interval ]Φ where H(s) = 1 for s ≥ 0 and H(s) = 0 otherwise, J cr,α is a critical value for J α , while U 0α , q and r are material parameters [95]. We remark that, in principle, q, r are different in the aggregate and in the plates. Nevertheless, in this work, for the sake of simplicity, they are taken to be equal in the aggregate and in the plates. Since in the plates, i.e., for β = u or β = , the solid phase is much stiffer than that of the aggregate, we prescribe a Helmholtz free energy density of De Saint-Venant type, i.e., where λ β and μ β are Lamé's constants of the βth plate and, since no remodelling occurs in the plates, coincides with the elastic Green-Lagrange deformation tensor. Hence, by accounting also for (33), the constitutive expression of the second Piola-Kirchhoff stress tensor of the solid phase of the βth plate is given by For completeness, we recall that the constitutive parts of the first Piola-Kirchhoff stress tensor and of Cauchy stress tensor are related to each other through P β , where, with a slight abuse of notation, σ (sc) β is here understood as a function of material points and time. We notice that the use of the penalty potential U β for the plates (β = u, ) is not mandatory as long as the plates are stiff enough to avert compaction. Nevertheless, there can be situations in which the stiffness of the plates becomes comparable with that of the aggregate, as is the case, for instance, when the contact between tissues is simulated. For this reason, and since the implementation of the penalty potentials constitutes only a minor computational caveat, we prefer to formulate a general model, even though, in this work, compaction does not arise and, thus, neither the penalty potentials nor their related stresses are activated.
Since the aggregate is assumed to be soft, isotropic and hyperelastic, we describe its mechanical response by means of a Holmes&Mow strain energy density function [53], i.e., where α 0 , α 1 , α 2 and α 3 are material parameters, such that it has to be α 3 = α 1 + 2α 2 , and I 1e , I 2e and I 3e are the three principal invariants of C e , i.e., so that the strain energy density can be rephrased as a function of C a and B p : Consequently, the constitutive part of the second Piola-Kirchhoff stress tensor of the solid phase in the aggregate is given by is expressed entirely in terms of C a and B p . Also in this case, the constitutive parts of the first Piola-Kirchhoff stress tensor and of Cauchy stress tensor satisfy the chain of equalities P (sc) is here a function of material points and time. Following [9,42,46,89], the plastic-like metric tensor B p can be assumed to vary according to a flow rule of Perzyna-type, given by one of the two equivalent forms (see also [9,46,89]) both of which are compatible with the part of the residual dissipation inequality associated with remodelling, which reads D rem = σ (sc) a : g p = Σ (sc) a : gL p [46] (see also the previous works of Armero, e.g. is the constitutive part of the Mandel stress tensor, L p =Ḟ p F −1 p is the rate of the plastic-like distortions and p = F e L p F −1 e is its elastic "push-forward". In Eqs. (40), γ p is given by [40,42,46], where λ p is a positive phenomenological parameter accounting for the time scale at which remodelling occurs, σ y is the yield stress of the aggregate and [s] + = 1 2 (|s| + s), for any real number s. Further details regarding the two equivalent forms of the remodelling law in Eqs. (40) can be found in [46]. Here, we simply mention that a more accurate definition of γ p should include the dependence of λ p and σ y on the a (x, t) = 0, no remodelling occurs at that point (see e.g. [40,41]). Unless stated otherwise, throughout this work we assume the spatial permeability tensor to be "unconditionally isotropic" (in the jargon of [11]), so that we can write and, with a slight abuse of notation, where the constitutive expression of k α is of the Holmes&Mow type [53], i.e., with m 0 and κ being material parameters, and k 0α being a scalar reference permeability [53]. We remark that, in principle, m 0 and κ may be different in the aggregate and in the plates. Nevertheless, in this work, for the sake of simplicity, they are taken to be equal. It is worth to remark that assuming the same deformation-dependent functional form of the scalar permeability k α both for the aggregate (α = a) and for the plates (α = u, ) may seem an unnecessary complication of the model. Indeed, the plates are much less deformable than the aggregate and would thus admit a permeability independent of deformation. However, there are at least two advantages that make us prefer the present formulation of the model. The first one is that we have a rather general and flexible approach, which can be easily adopted to different situations in which the plates are not so stiff as in the experiments that we have simulated.
The second advantage is that it is sufficient to implement a single functional form of the permeability and then associate it with different computational domains (i.e., the multicellular aggregate and the plates) by simply specifying the corresponding material parameters. On the other hand, the generality and flexibility of our formulation do not lead to an appreciable computational effort. Indeed, the relatively small deformations in the plates are such that the linearisation procedure required by the deformationdependent permeability employed in our work is not so computational onerous as it may seem. This is due to the fact that small deformations only need few iterations of the nonlinear Newton solver to reach convergence within each time step. Assuming the permeability tensor to be "unconditionally isotropic" [11], as in Eq. (42), implies that k α is spherical, which means that its principal axes are influenced neither by the deformation nor by the remodelling. Especially for this independence on remodelling, we shall be referring to the permeability tensors defined in Eqs. (42) and (43) as remodelling-independent permeability tensors. However, as anticipated in Introduction, in addition to Eq. (42), we also consider the case in which the permeability tensor of the aggregate is a constitutive isotropic tensor-valued function that does not reduce to a spherical tensor. Yet, we maintain the hypothesis of purely spherical permeability tensors for the plates.
By using the Representation Theorem as indicated in [11], and adapting the expression reported therein to our notation and, above all, to the characterisation of isotropy with respect to the aggregate's natural state, we specify k a also through the constitutive law where k (1) a and k (2) a are scalar-valued constitutive functions given by (cf. [11]) with k (1) 0a and k (2) 0a being material parameters regarded as constants in the present framework. More specifically, in the sequel, we prescribe k (1) 0a and k (2) 0a to be where c is a strictly positive constant that will be determined by means of a best-fit procedure. The purpose of introducing the constitutive law (45) is to study its consequences on some of the most relevant results of our model, as will be discussed in Sect. 7.4.
To emphasise how the deformation gradient tensor and the plastic-like distortions contribute to the constitutive expression of the permeability tensor defined in Eq. (45), we consider the Piola transformation of k a , i.e., Because of the influence of remodelling on the permeability tensors defined in Eqs. (45) and (48), we shall denominate them remodelling-dependent permeability tensors.
All the material parameters introduced up to here are collected in Table 1.

Summary of the model equations and benchmark test
To investigate the capability of the aggregate of recovering its original shape in spite of its structural changes and of the other dissipative processes occurring in it, we simulate an experimentally relevant uni-axial compression-release test. We keep the upper plate fixed over the whole time interval T , so that the aggregate is compressed by moving the lower plate. More specifically, the aggregate is compressed according to a specific loading history and it is then unloaded so as to observe how it relaxes (e.g. how the stress evolves in time) and to which extent it recovers the shape it had prior to the loading. We summarise the model equations in Eulerian form, which read We assume the contact between the aggregate and the plates to occur in the absence of friction. Hence, we do not prescribe any condition along the tangential direction of the contact areas Σ u (t) and Σ (t), and we impose the following conditions [52] with n β being the field of unit vectors normal to Σ β (t) and β ∈ {u, }. In particular, Eqs. Σ u (t) and Σ (t) and of the normal component of the solid velocity. Moreover, Eq. (50d) expresses the continuity of the mass flux across Σ β (t). In addition, we consider the boundary conditions where ∂ T B u (t) denotes the top surface of the upper plate and ∂ B B (t) the bottom surface of the lower plate. For future use, we introduce also the notation ∂ T B (t) to indicate the top surface of the lower plate. Note that, for the sake of a lighter notation, in Eqs. (50a)-(50d) and (51a)-(51f), the dependence on spatial points and time is omitted. In Eqs. (51b)-(51f), n α is a field of unit vectors normal to the boundary of B α (t), with α ∈ {u, a, }, except for a set of zero measure in which the unit normal vector is undefined. Moreover, Eq. (51a) describes the fact that the fluid phase is in hydrostatic equilibrium with the fluid of the bath in which the aggregate is immersed, Eqs. (51b)-(51d) are null traction conditions, whereas Eqs. (51e) and (51f) express that the boundaries Γ u (t) and Γ (t) are impermeable. Finally, we supply Dirichlet conditions on the motion in Lagrangian form as Equation (52a) is a zero displacement condition, whereas Eq. (52b) provides the expression of the imposed displacement u(t), with e z being the unit vector aligned with the z-axis. In particular, following [42], the loading history is prescribed as where ]0, t ramp [ is the time interval over which the lower plate is moved upwards to compress the aggregate against the upper plate, which is kept fixed; [t ramp , t end − t ramp [ is the time window over which the load is kept constant (the aggregate remains compressed); [t end − t ramp , t end [ is the interval over which the load is gradually removed; [t end , t end + t rest ] is the time window over which the relaxation and the shape recovery are evaluated. Note that we exclude t = 0 s from (53) because the configuration of the multicellular aggregate prior to the launching of the simulated experiment (i.e., before triggering the displacement of the lower plate) is not the undeformed reference configuration of the aggregate. Rather, it is a configuration that is already deformed because of the force of gravity until time t = 0 − s, to which we are superimposing the prescribed displacement from t = 0 + s on. In passing, we notice that the final time of the simulation, i.e., the instant T introduced in the definition of T = [0, T ], is given by t end + t rest . By denoting by ∂∂ T B u (t) and ∂∂ T B (t) the circumferences representing the boundary of ∂ T B u (t) and ∂ T B (t), respectively, the boundary conditions are completed by imposing, for the pressure, Dirichlet points on ∂∂ T B u (t) and ∂∂ T B (t). On these circumferences, the pressure is set equal to zero. This stabilises the simulations.
To conclude the formulation of the problem, we mention that the initial conditions are given in Lagrangian form, along with the equations in the bulk, i.e., which are the ones that are solved numerically. We recall that, in the reference configuration, the system is undeformed, with zero pressure, and finds itself in a state characterised by the absence of plastic-like distortions. Hence, the initial conditions read The computational issues related to the numerical implementation of the contact conditions will be discussed in the following section.
The parameters that describe the physical and geometrical properties of the multicellular aggregate are reported in Table 1. In particular, since the simulations are run for different values of the aggregate's shear modulus μ a , and Poisson's ratio ν a , we express the material parameters α 0 , α 1 and α 2 as functions of μ a and ν a , i.e., where α 3 is set equal to unity [95].
We emphasise that our model, which describes the aggregate as an elasto-plastic porous medium saturated by a fluid in Darcian regime, naturally introduces two independent characteristic time-scales: one is dictated by the permeability in Eq. (49a) and the other one is due to the evolution law of the plastic-like distortions in Eq. (49c). We remark that also the paper by Forgacs et al. [32] features two characteristic times, determined by solving a second-order ordinary differential equation in the force exchanged between the aggregate and the plates, and leading to a bi-exponential law of the relaxation of such force. However, the modelling approach followed in our work and, more specifically, the physics that we are describing are deeply different from those in [32]. This fact, in turn, is reflected, among other aspects, in the different nature of the characteristic time-scales introduced by our model.
The contact problem was addressed by defining two Contact Pairs, one for I u and the other one for I , and adding the contact feature with the most appropriate boundary conditions, among those available in COMSOL TM (e.g. boundary conditions of Dirichlet or Neumann type), in the packages Structural Mechanics and Darcy's Law [18,19]. On the other hand, we did not select any contact specific option for the ODE associated with B p because no boundary conditions are needed for it [18,19].
To ensure the correct resolution of the contact between the aggregate and each plate of the apparatus, which constitutes the most peculiar aspect of the simulations, we had to use some precautions in our  set-up. In fact, to make the simulations start, we had to introduce a "critical" distance between each plate and the aggregate, with the purpose of "telling" the software the threshold below which the contact was to be detected. This was done in the package Structural Mechanics and, for coherence with the explanation given in the COMSOL TM manual [19], we chose the value h = 10 −1 μm for the distance between the contact surface and the reference surface. Note that, since this choice yields h/R = 10 −3 , it guarantees a "good" scale separation with respect to the radius of the aggregate R = 100 μm (see Table 1) and, thus, an acceptable "advance" of the contact dynamics. In our simulations, an Augmented Lagrangian method for the solution of contact has been employed [42,90].
To reduce the numerical errors, the mesh sizes for the plates and for the aggregate were chosen by setting the options "Extremely Fine" and "Extra Fine", respectively, which means that the mesh in the plates is finer than that in the aggregate. Indeed, these mesh sizes prevent the onset of abnormal pressure oscillations in the contact areas and make sure that the pressure in the neighbourhood of such areas does not acquire unphysical values [18,19]. For example, for different mesh sizes we observed values of pressure in the proximity of the contact areas that were orders of magnitude higher than the values of the pressure in the contact areas. Moreover, according to the COMSOL TM manual [19], the mesh size used for discretising the plates is approximately half the mesh size used for the aggregate. Finally, after several tests, the maximum time step was set equal to Δt = 0.125 s, because higher values gave rise to numerical errors in the resolution of Eq. (54a) in the plates.
We remark that, although the geometry of the problem at hand suggests that contact should start in the points of the aggregate and of the plates lying on the symmetry axis of the system as a whole, in our simulations COMSOL TM seemed to be unable to catch this feature automatically. Hence, to circumvent this difficulty, we proceeded in two steps. First, we called for the contact control variables ''incontact p1'' and ''incontact p2'' [18,19], which represent the indicator functions of the contact surfaces, and whose indices ''p1'' and ''p2'' refer to the first and to the second pair of the contact pairs created in our model. Then, we modified ''incontact p1'' and ''incontact p2'' by requiring that the points of the boundary of the aggregate and, thus, the associated points of the plates in the cylindrical neighbourhood of the axis of symmetry with radius of r = 10 −3 μm remain in contact until the compression is relaxed. Again, note that the ratio r/R = 10 −5 ensures that the alteration of ''incontact p1'' and ''incontact p2'' affects mesh points lying so close to the system's symmetry axis that the simulations are not appreciably influenced by this additional contact condition.
Note that COMSOL TM features an option to scale automatically the unknowns that have to be computed. In the case in which remodelling is active, the automatic scaling supplied by COMSOL TM is sufficient for a correct resolution of the pressure field. On the contrary, in the absence of remodelling, the scaling provided by COMSOL TM was not sufficient and it was, thus, necessary to conduct a manual scaling. Hence, in addition to the previous precautions, for the simulations in which the equation for B p (54c) is deactivated (no remodelling, but unchanged boundary conditions, parameters and geometry), a manual scaling for the pressure p of 10 MPa is required. These simulations are essential to compare the different phenomena that are triggered by the plastic-like distortions.

Results and discussion
In this section, we discuss the results of the above-described compression-release test that are most relevant to us and we propose an interpretation of the related physical phenomena occurring in the aggregate. As anticipated above, the main purpose of this work is the achievement of a deeper understanding of the physics concealed in the presented model through a thorough visualisation of the system's evolution. In particular, we will focus on the description of the behaviour of the fluid phase and on the interactions between the evolution of the pressure distribution and the plastic-like distortions.
Before discussing the results, we mention that, in Sects. 7.1, 7.2 and 7.3, the spatial permeability tensor of the aggregate is taken equal to the spherical tensor defined in Eq. (42). We do this in order to highlight how the interplay between the evolution of the fluid and the remodelling of the aggregate is visible, although less direct, also in the case of a remodelling-independent permeability tensor. In Sect. 7.4, instead, we will study separately the case of remodelling-independent permeability from that of remodelling-dependent permeability. This is done to emphasise how a more direct effect of the remodelling on the fluid motion may influence the tendency of the aggregate to restore its original shape after compression.

Pressure
Since our model of fluid motion falls within Darcy's regime, the fluid velocity can be reconstructed from the pressure distribution as predicted by Darcy's Law, i.e., v However, the role of gravity, represented by (f) α ga g , turns out to be negligible in comparison with the pressure gradients generated by the compression. For this reason, the fluid motion is essentially driven by the pressure gradients alone. To emphasise this, we notice that the initial pressure distribution is not uniform because of gravity, although its evolution will not be significantly influenced by (f) α ga g .

The role of contact on the pressure distribution.
In the initial configuration of the system, i.e., before the prescribed displacement starts, the aggregate, which is placed on the lower plate, is stressed because of gravity. This situation generates an overpressure in the proximity of the contact zone of the aggregate with the lower plate, as a consequence of the weight of the aggregate itself.
During the application of the prescribed displacement, the lower plate moves upwards and compresses the aggregate against the upper plate, which remains fixed, and an overview of the spatial distribution of pressure at some specific instants of the loading is reported in Fig. 3. In particular, in correspondence to the regions near the plates, the pressure features a spike, whose amplitude grows gradually in time (coherently with the speed of the loading ramp) as the aggregate's porosity diminishes (see Fig. 4 for t = t ramp = 5 s). This is accompanied by an increase of the contact area of the aggregate with the plates.
The behaviour just described is largely influenced by the way in which the contact progresses and by the permeability of the plates. To see this, let us recall that the outflow of the fluid through the free boundary of the aggregate is entirely due to the aggregate's permeability and to the pressure gradient, whereas the passage to the plates depends also on their permeability. Given this situation, if the permeability of the plates were comparable with that of the aggregate, or even higher, the fluid would migrate easily from the aggregate to the plates, and its pressure would remain relatively low. Yet, since the plates are much less permeable than the aggregate, they limit the transit of the fluid, thereby triggering an increase in the fluid pressure in the aggregate.
We remark that the assumption of a relatively low permeability of the plates modelled in our simulations is meant to render the simulated experiment as close as possible to the real one, in which the plates are impermeable (see Sect. 2.2). Clearly, in this latter case no passage of fluid from the aggregate to the plates would be observed.

The "syringe effect".
As the displacement of the lower plate is maintained, the fluid follows the pressure gradients and it distributes in the central area of the aggregate, around the vertical axis of symmetry, as can be seen in Fig. 3b. During the time interval ]0, 40[ s, which comprehends the application of the prescribed displacement, a loss of mass is registered, as an outflow of fluid occurs through the free boundary of the aggregate in response to compression (see Fig. 3a-c). This can be also deduced by looking at the isobaric lines drawn in Fig. 3. Indeed, with reference for instance to Fig. 3b, c, one observes that the isobaric lines descend from a high pressure core towards the free boundary over which the pressure is atmospheric. Therefore, coherently with this pressure distribution, and according to Darcy's law in the case of a spherical permeability tensor and negligibly small role of gravity, the velocity of the fluid relative to the solid is parallel to the pressure gradient and, thus, orthogonal to the isobaric lines. Before proceeding, we notice that, in the first instants of the interval ]0, 40[ s, close to the boundary (i.e., at X U = (40,0,190) μm of the reference configuration as indicated by the black spot in Fig. 3), and regardless of whether or not remodelling starts, the fluid pressure increases and the porosity decreases (see Fig. 4).
The behaviour of the aggregate depends noticeably on the activation of remodelling and on the concomitant redistribution of the stresses and of the fluid pressure. In particular, the evolution of the latter one triggers an effect that is not visible in the absence of remodelling and that, even in the presence of remodelling, emerges only below a certain threshold initial diameter of the aggregate. After the initial increase of pressure and decrease of porosity, this effect manifests itself through two aspects: (i) a fall of the fluid pressure below zero and its subsequent ascent towards zero from 40 s onwards; (ii) a growth of porosity over the same time window. As a consequence of these two facts, the fluid tends to flow back into the aggregate, since it finds enough space due to the reopening of the pores (see Fig. 4). This behaviour was termed "syringe effect" in [21,23].
We emphasise that the "syringe effect" was observed in the completely different situations studied in [21,23], where it was ascribed to the combination of remodelling and anisotropy in fibre-reinforced anisotropic media. Indeed, no such effect could be reckoned in the isotropic cases addressed in [21,23], even in the presence of remodelling. Moreover, a similar effect was also observed in [81] in concomitance to the introduction of non-local effects in the context of tumour growth. In our setting, instead, the "syringe effect" occurs in the isotropic aggregate, only in the presence of remodelling, but exclusively below a given characteristic length of the diameter of the aggregate prior to the compression-release test. Probably, this is due to the fact that, over a certain aggregate size, the backflow of the fluid is dispersed in a bigger volume, so that the "syringe effect" is no longer appreciable. This suggests us to interpret the "syringe effect" as a possible size effect, at least in the present context, rather than as an effect directly related to anisotropy, which, in turn, could magnify the effect itself. This requires further investigations.

The "pressure bubble".
Because of the plastic-like distortions (see Fig. 5), the evolution of the pressure manifests itself as a travelling "pressure bubble" that accompanies the "syringe effect", thereby reflecting the inversion of direction of the fluid flow in the aggregate. Indeed, for t ≥ 30 s, the interstitial fluid migrates towards the contact areas, which will be then characterised by negative pressures. Such "pressure bubble" is a subset of B a (t) in which the fluid pressure is positive. We remark that the formation of the "pressure bubble" is a characteristic feature of remodelling (see Fig. 3), since no reopening of the pores and no formation of negative pressure areas occur in the purely elastic biphasic model. In other words, there is a tangible difference between the two considered models, i.e., the one with and the one without the account for structural transformations (remodelling): in the first case (with remodelling switched on), after the "pressure bubble" disperses while approaching the free boundary of the aggregate Γ a (t), no fluid leakage is observed on Γ a (t) and, rather, an inflow of fluid is registered through this surface; in the second case (no remodelling), the "pressure bubble" does not arise and the fluid continues to leak through Γ a (t).
It should be noted, however, that, even in the presence of remodelling, the "pressure bubble" does not form when the elastic parameters of the solid phase are relatively small, as is the case for μ a = 3 kPa, ν a = 0.2. This is because, for such values of the elastic parameters, the constitutive part of the mechanical stress is only slightly greater than 2/3 σ y and, thus, remodelling is barely activated and evolves slowly because γ p is rather low.
Coherently with these comments, we remark that the "pressure bubble" evidences the tendency of the remodelling to lower the magnitude of the stresses, and thus also of the fluid pressure [23,46] (which, in fact, becomes negative around the bubble) in response to the applied load. On the other hand, the nature of such stress relaxation is very different from the one originating in the aggregate during the  Table 1. The isobaric lines are perpendicular to the portions of the boundaries of the upper and lower plates that are not in contact, coherently with the no-flux conditions. The "pressure bubble" visibly forms only for t > 35 s in the case with remodelling unloading phase, which occurs regardless of whether its remodelling is set off or not. This brings us to another consideration: after the unloading phase, the aggregate is characterised by negative pressure independently of remodelling, which means that the fluid is allowed to flow back into the aggregate. Also this result contributes to shed light on the dynamics of the interstitial fluid, since it underlines that, in the presence of remodelling, the fluid enters the aggregate even during the loading phase, thereby giving rise to a behaviour that is not observed in the purely elastic case, as shown in Fig. 3.

Remodelling
Following [23,42,46], also in this work the onset and evolution of plastic-like distortions are stress-driven processes that take place when the Von Mises stress of the solid phase, σ VM = 3/2 ||devσ (sc) a ||, exceeds the yield stress of the aggregate, σ y . In all the other cases, the solid phase of the aggregate behaves elastically. In fact, depending on the size of the aggregate and of the model parameters, the plastic-like distortions may exist already in the initial configuration B a (0), if the force of gravity (i.e., the only load prior to the compression) induces a sufficiently high Von Mises stress.
In our compression-release test, the commencement of remodelling is observed already during the initial loading: σ VM spikes at t = 5 s, but it diminishes as the change of the internal structure, triggered by the plastic-like distortions, makes the aggregate relax (Fig. 6). Moreover, the presence of the fluid contributes to lower the internal stresses, but the pressure does not interact directly with B p , since only  (40,0,190) μm in the reference configuration of the aggregate of initial radius R = 100 μm (curves marked with the circles and the triangles). In the case of remodelling and radius R = 100 μm (line marked with circles), the pressure equals zero at t = 40 s. We consider the constant ratioũmax := umax/2R = 0.3, so that umax = 60 μm. For an aggregate of initial radius R = 900 μm, same constant ratio and, thus, umax = 540 μm, we select X U = (360, 0, 1710) μm (curve marked with the squares). The other parameters are reported in Table 1 the deviatoric part of the constitutive Cauchy stress tensor, i.e., σ (sc) a , features in Eqs. (40). In addition, we notice that, during the compression phase, the Von Mises stress attains its highest value in the centre of the aggregate and the lowest one in the contact zones. However, it inverts its distribution after the plates are moved away (see Fig. 6). For t = 40 s, the regions of the aggregate that were in the proximity of the contact boundary feature the highest value of the Von Mises stress, while, in the central region, the Von Mises stress is one order of magnitude lower. It should be noticed, however, that the magnitude of the Von Mises stress near the contact surfaces passes from a value close to 4.0 × 10 3 Pa before the detachment of the aggregate from the upper plate to a value of about 3.5 × 10 3 Pa attained after such detachment. This inversion is supported by two aspects of the model: the presence of plastic-like distortions and the spatial distribution of porosity. As visible in Fig. 5, almost no plastic-like distortions take place in the proximity of the contact boundaries and, thus, the effect of the relaxation of the internal structure is much lower. In addition, the pores near to the contact boundaries "feel" the compression more intensively (Fig.  4) than in the central region. Indeed, during the removal of the compression, the central region unloads faster (in the time window ]35, 37[ s) because, although it features higher values of the Von Mises stress, the pores are less compressed. So, for t = 40 s, the regions that were compressed are still under the elastic response and the distribution of the Von Mises stress is reverted with respect to the compression stages (see the transition between Fig. 6a, d).
With reference to [42], in which the compression-release test is simulated with the same loading history, we notice that the distribution of the Von Mises stress after the removal of the compression is recovered in our work for t = 60 s. In fact, this is an expected behaviour as the introduction of the fluid phase increases the time needed for the shape recovery with respect to the monophasic model, so that, for t = 40 s, the biphasic model is still undergoing a notable expansion. For t = 60 s, the shape recovery is almost completed and a similar distribution of the Von Mises stress is recovered.
Finally, the evolution of the stresses in the aggregate is accompanied by the evolution of the plasticlike distortions, which we describe here through the Frobenius norm of the plastic Green-Lagrange strain tensor E p := 1 2 (C p − G), i.e., E p = (E p : E p ) 1/2 . It is important to notice that, as the simulation progresses, the magnitude of the Von Mises stress diminishes (see Fig. 6) and E p evolves towards a stationary distribution (see Fig. 5). In fact, the plastic-like distortions, after their onset, accumulate during the compression stages and tend to remain permanent afterwards. We also note that, for t = 60 s,   [42]. In these pictures, we took tramp = 5 s and t end = 40 s. The difference in the distribution of the Von Mises stress between the plates is due to the presence of gravity i.e., for a time instant in the interval [t end , t end + t rest ], the stress distribution is still close to that of the material in the relaxation phase, i.e., that of the phase subsequent to the end of the maintenance of the compressive displacement, i.e., t = t end − t ramp = 35 s. Also in this case, the influence of the loading history on the mechanical state of the aggregate is measured by the plastic-like distortions, which, as already mentioned, tend to a stationary distribution. Therefore, if the stationary configuration of the aggregate were taken as the new "initial" configuration for another experiment, one should consider an initial configuration with residual distortions.

Stress relaxation curves
In this subsection, we present a qualitative comparison between the experimental stress relaxation curves taken from [8,32,55] and those obtained in our simulations for the fully biphasic model. To calibrate our model, we vary the parameters in Table 1 over a range of values compatible with those reported in the literature, and we are able to retrieve qualitatively the target experimental curves for physically acceptable values of the parameters. This, in turn, permits us to appreciate the stability of the model within the chosen scales of the parameters. Before going further, we emphasise that we did not conduct a direct quantitative validation of our numerical predictions with the experimental curves in [8,32,33,55], since this would require further biological data and details on the mechanical tests that are not available in [32]. However, we showed that our model is able to reproduce a sufficiently wide range of normalised stresses by varying the model parameters, so that different cellular populations can be described by changing the combination of the parameters of our model, possibly supported by other mechanical tests.
We consider the stress relaxation curves obtained by varying the remodelling parameters σ y and λ p , the elastic parameters μ a and ν a , the scalar permeability of the aggregate k 0a , the maximum normalised displacementũ max , and the initial radius of the aggregate.
The stress relaxation curves presented in Fig. 7 are obtained by integrating the negative of the normal component of the total Cauchy stress tensor σ a = −p a g −1 + σ (sc) a over the contact area Σ u (t). We call the resulting force total exchanged force. On the other hand, the stress relaxation curves in Fig.  8 are obtained by integrating only the negative of the normal component of σ (sc) a . Accordingly, we conventionally denominate the resulting force constitutive exchanged force. A comparison of the total exchanged force with the constitutive one is helpful to enlighten the overall influence of pressure (see Fig.  9). In Figs. 7, 8 and 9, the curves represent how the exchange of forces through the contact area Σ u (t) varies in time. Note that the terminology "normalised force" reported in these figures refers to the fact that the force (expressed in mN) is divided by the value of the gravity acceleration a g = 9.81 m s −2 , thereby yielding the units of mg. This is done to help the Readers in the qualitative comparison of our results with the experimental curves obtained by [32,33], in which the values of the force are given in mg.
With reference to Figs. 7 and 8, we notice that, in the presence of remodelling, after the initial peak corresponding to the applied displacement, the onset of plastic-like distortions causes a progressive reduction in the exchanged force until a stationary value is asymptotically achieved. In the case in which remodelling is not active, so that the behaviour is purely elastic, a stationary value is found shortly after the application of the prescribed displacement (data not shown for the purely elastic case).
As reported in Fig. 7a, the yield stress is varied while maintaining an imposed displacement ofũ max = 0.3. This value agrees with the experimental observation of Jakab et al. [55], who observed that aggregates may be compressed up to a maximum of "approximately 30% of their original diameter" in order to prevent cells from the occurrence of damage and intense shape modification of the cellular structure. We remark that a change in the yield stress does not vary significantly the amplitude of the initial peak of the total exchanged force shown in Fig. 7a. On the contrary, the stationary value of the total exchanged force is influenced by the value of the yield stress, which is in agreement with the results obtained in the one dimensional analysis reported in [40]. Indeed, a lower value of σ y triggers the remodelling for lower values of the Von Mises stress, thereby facilitating the commencement and the accumulation of the plastic-like distortions. Also, the amount of exchanged force due to σ (sc) a exhibits a similar behaviour, although the initial peaks are positioned slightly forward and the stationary values in Figs. 7 and 8 coincide. In addition, for the considered set of parameters and for the type of physical interactions accounted for, the force exchanged through the contact surface depends mostly on the constitutive contributions of Cauchy stress. On the other hand, the influence of the pressure is essentially restricted to the initial peaks, as is clear from Fig. 9a, where the force exerted by the pressure over the plates becomes very small after relatively short time intervals. We remark that the influence of the pressure is in percentage more relevant for larger prescribed displacements (see Figs. 7e and 9e). Indeed, the amplitudes of the peaks are higher in the relaxation curves of the total exchanged force (Fig. 7e) than in the ones obtained for the constitutive exchanged force (Fig. 8e).
The remodelling parameter λ p and the yield stress σ y are hardly detectable because there are various factors that can influence their experimental determination. In this sense, the comparison of the simulated stress relaxation curves with the experimental ones may be helpful to identify the values of λ p and σ y that best fit the experimental data, when such a procedure is done. As shown in Fig. 7b, varying λ p plays a relevant role in the determination of the height of the peak and of the shape of the relaxation curves, since λ p is related to the characteristic time of remodelling [21,41], whereas it does not affect appreciably the stationary value of the exchanged force. Note that, for the curves corresponding to the values λ p = 0.1 (MPa s) −1 and λ p = 0.5 (MPa s) −1 the stationary values are not reached in Fig. 8b. In particular, if we compare the relaxation curves in Fig. 7b with the ones in Fig. 8b, it appears that λ p has a  Table 1. We emphasise that all the curves plotted in this figure refer to the case with remodelling role on how the pressure varies in response to remodelling (Fig. 9b). For increasing λ p , remodelling occurs faster and, accordingly, the stress magnitude diminishes faster too: for fixed time, the reduction of the stress magnitude is more marked. This is more visible for the relaxation curves associated with σ (sc) a (see Fig. 8b), whereas it is less visible for the curves associated with σ a = −p a g −1 +σ (sc) a (see Fig. 7b) because, with respect to λ p , the pressure peaks decrease slightly more slowly than those associated with σ (sc) a . Moreover, looking at Fig. 8b, we notice that the difference of the peak attained for λ p = 1.0 (MPa s) −1 and the one attained for λ p = 0.1 (MPa s) −1 is about 4 mg, while it is about 2 mg in Fig. 7b. In saying this, we have accounted for the fact that the two peaks are attained at different instants of time.
As done in [42], we study also the sensitivity of the model to the variation of the elastic parameters by considering a biologically admissible range [76,85,94]. In this respect, we notice that, if Poisson's ratio is kept fixed and the shear modulus is increased, or vice versa, the magnitude of the initial peak of the total exchanged force increases too (Fig. 7c). The same occurs also for the constitutive exchanged force (Fig. 8c) and pressure (Fig. 9c). All these results are consistent with the monophasic model presented in [42]. To us, this is an indication that our modelling framework and simulations are reasonable, since the presence of the fluid phase changed significantly the model equations and the boundary conditions adopted in this work as well as the computational effort for their solutions. The variations of the mechanical parameters of the model allows to take into account the high variability of the stress measurements reported in the literature for different cell lines [32,33] and for slightly different experimental settings [8,55]. For instance, the curves  Fig. 7b marked with circles and squares). The amount of stress exchanged between the aggregate and the plates was observed to vary nearly sevenfold in [55], depending on the amount of latrunculin concentration, that affects the adhesive and mechanical properties of the cells inside the aggregate. We remark that such variations in the values of cell mechanical features are not unusual in the biological works: the mechanical parameters are strongly dependent on the cell type considered and, even considering the same cell line, the mechanical properties may range over several orders of magnitude in response to the time-scale of the measurements, applied stress and deformation rate [76,101].
With the introduction of the fluid phase, and because of the hypothesis that the plates are saturated porous media, a permeability tensor has to be introduced both for the aggregate and for the plates. As anticipated in Sect. 2.2, although at first sight this approach may seem to lead to superfluous complications, it is supported both by technical reasons and by the flexibility of the mathematical model. In fact, its advantages from the point of view of modelling have been already discussed in Sect. 4.4, and we thus focus here on the technical motivations at the basis of our framework. In this respect, we report that, in some  Table 1. Note that, also in this case, all the curves plotted in this figure refer to the case with remodelling simulations in which the plates of the experimental apparatus were taken as impermeable monophasic solids (data not shown), we faced a problem that is inherent in the occurrence of contact. This problem manifested itself through the difficulty of simulating the contact between the impermeable solid plates and the aggregate, which in our model is a saturated porous medium, i.e. a fluid-solid mixture at the considered scale of observation. To our knowledge, the problem of the contact between a mixture of this type and a monophasic solid is not straightforward already at the stage at which the correct boundary condition for the normal component of the stress tensor has to be prescribed on the contact surface. The complication, indeed, arises because the total stress in the mixture features a component due to the fluid pressure, which is absent in the monophasic solid (see e.g. [52]). In our study, this description becomes even more involved because of the presence of remodelling, which complicates the equations and affects the stress distribution in the aggregate. However, beyond these physical complications, in the default settings of COMSOL Multiphysics TM v. 5.3a, the implementation of contact seems to require to define the same variables in the two domains that enter in contact with each other, because each domain accounts solely for the variables that are defined within it. Then, when the COMSOL TM command Contact Pair is invoked, it connects reciprocally the variables that describe the same physical entity, associated with the different domains of the considered contact pair. In our case, the software would initialise the pressure in the aggregate, but, if the plates were not modelled as saturated porous media too, the software would not "see" an allocated variable associated with the pressure field in them. To circumvent this technical obstacle in a way that can be proficuous also from the side of modelling, we have assumed that the plates are fluid saturated porous media, and, to approximate the real experimental set-up, we have taken their permeability at least four orders of magnitude smaller than that of the aggregate. It is important, then, to keep this difference in the orders of magnitude fixed, even when the permeability of the aggregate is varied to assess its repercussions on the mechanic and hydraulic behaviour of the overall system. In particular, we notice that, for decreasing permeability, the peaks of the total force exchanged through the contact surface become higher (see Fig. 7d). This is because the fluid exerts a higher pressure on the contact surface. Still, at the same time, the force associated with the constitutive contribution exchanged through the contact surface decreases significantly (see Fig. 8d). Conversely, both the characteristic relaxation time and the stationary value of the total exchanged force are almost unaffected by permeability changes. We notice that a behaviour similar to the one shown in Figs. 7d and 8d is also observed in Fig.  9d.
Beyond the parametric study of the aggregate's mechanical response for varying remodelling parameters σ y and λ p , elastic coefficients μ a and ν a , and reference permeability k 0a , it is also important to assess how the aggregate's behaviour is influenced by the imposed external stimuli. Since in the present framework the external stimuli are supplied by a prescribed displacement, whose magnitude must be related to the initial size of the aggregate through its radius, we study the parametric response of the aggregate by varying the ratioũ max = u max /(2R) (Figs. 7e, 8e and 9e).
Lastly, a variation of the geometric parameters is considered and, in agreement with [42], it is observed that increasing the radius of the aggregate intensifies significantly the contact force both at the initial time and when the stationary state is approached (Figs. 7f, 8f and 9f). We notice, in particular, that the pressure contribution becomes more relevant than that of the constitutive exchanged force for increasing aggregate's radius (see Fig. 9f). Similarly, as the prescribed displacement increases, the contribution of the pressure becomes more relevant, until it reaches a point after which it exceeds the constitutive one.
We remark that the sensitivity analysis shows that the radius of the aggregate and the imposed deformation highly influenced the amount of the stress exchanged between the aggregate and the plates. Therefore, in order to qualitatively compare the numerical results with the experimental ones, such values should be precisely known (whilst they are not always reported in the literature).
We highlight that in Fig. 9a-f it is also possible to appreciate the "syringe effect", which manifests itself as a change of sign of the force related to the pressure exchanged through Σ u (t). Indeed, while the lower plate is pushing upwards, so that the aggregate is compressed and, thus, contracts, the inflow of fluid due to the "syringe effect" tends to inflate the aggregate, which, consequently, pushes on the plates (specifically, in Fig. 9a-f, we observe the effect of this pushing on the upper plate).

Shape recovery curves
While the shape recovery for compressed multicellular aggregates has been studied in [42] for the monophasic case, in our work we study whether the interstitial fluid has a role on the tendency of the examined aggregate to restore its original shape after compression. More specifically, we investigate whether the irreversible changes of the aggregate's internal structure are modulated by the fact that the fluid absorbs part of the stress induced by the applied load. To this end, since these irreversible changes prevent the aggregate from recovering its original shape, we observe the difference between the final and the initial shape of the aggregate and we relate it to the presence of the fluid through the permeability.
Both in the experimental procedures and in our simulations, measurements of the height and of the width of the aggregate are performed after a period of relaxation, during which the aggregate should recover its original shape if no plastic-like distortions took place. In fact, even in the presence of plasticlike distortions, a partial regain of the initial shape is observed, and the rate at which it occurs is slower than in the monophasic model and continues for a longer time. In our opinion, this is because the shape recovery is a two-level phenomenon for a biphasic aggregate: it involves, on the one hand, the reopening of the pores after compression and the return motion of the fluid, which flows back into the aggregate, and, on the other hand, the elastic part of the overall deformation, which is recovered after the release  Fig. 10. Each data point in the shape recovery curves represents the height/width ratio correlated with the time for which the aggregate is compressed. The first data point is for a compression time of 5 s. In the curve at the right the value of λp is set to be 0.5 (MPa s) −1 , since it is the value of the remodelling parameter for which the shape recovery curve is more similar to the experimental data extrapolated from [32] of the load. We remark that this behaviour is described in Fig. 6, as the regions proximal to the contact surfaces tend to recover elastically their original shape. In fact, the presence of the fluid allows to observe the shape recovery as though it were in slow motion. In the sequel, we consider separately the case of remodelling-independent permeability and that of remodelling-dependent permeability. Fig. 10, the final shape of the aggregate depends significantly on the value of the remodelling parameter λ p , since a faster onset of plastic-like distortions means that less compression time is required in order to affect the internal structure. The value λ p = 0.5 (MPa s) −1 is coherent with [42]. We remark, however, that also in this case the comparison is qualitative, since not all the data required to quantitatively compare the numerical results with the experimental ones are reported in the original paper by Forgacs et al. [32]. In fact, the qualitative agreement between our results and the experimental curve reported in [32] is good for short and long times of compression, i.e. when the experimental data are monotonically decreasing. Yet, for intermediate times, the experimental results are not well caught by the numerical predictions. Indeed, the experimental behaviour for intermediate times does not decrease monotonically, which, in our opinion, deserves further investigations in order to understand whether this behaviour is related to uncertainties in the measurements or to real biological features of the system. Lastly, the effect of the permeability variations, which is a dominant feature of the biphasic model, is more visible for periods of compression up to 5 min (see Fig. 10), with the value of the remodelling parameter kept fixed at λ p = 0.5 (MPa s) −1 . We notice that, for lower values of the permeability, the final shape of the aggregate is "less spherical". On the contrary, for higher values of the permeability, the aggregate tends to recover more efficiently its initial spherical shape because the interstitial fluid, escaped due to the compression, tends to flow back more easily into the aggregate. We also remark that, for longer compression times (i.e. for t end > 5 min), the fluid absorbs part of the stress and the shape recovery curve for k 0a = 2.5779 × 10 −4 mm 4 /(N s) returns to be close to the other ones.

Remodelling-dependent permeability.
In this section, we report on the influence of the use of Eq. (45) on the three major results of our work. To this end, we recall that the constant parameter c, by which k 0a is rescaled to obtained the permeability coefficients k (1) 0a and k (2) 0a (see Eq. (47)) is determined by means of a best-fit procedure of the shape-recovery curves. Fig. 11. Comparison among the experimental results extrapolated from [32], the numerical results obtained with the "remodelling-independent permeability" and the numerical results obtained with the "remodelling-dependent permeability". The wording "remodelling-independent permeability" refers to the spherical permeability tensor in Eq. (42), whereas "remodelling-dependent permeability" refers to the permeability tensor in Eq. (45). In the simulations, we set k 0a = 2.5779 × 10 −3 mm 4 /(Ns), k In particular, if, for example, we take the recovery curve in [32] for reference, then, by selecting a relatively short time window, we can find the value of c such that the computed height/width ratio best approaches the experimental data taken from [32]. We emphasise that, since not all the parameters used in our model coincide with those reported in [32], the determination of c is not meant to be quantitative. Rather, it is only meant to show how the introduction of a physically grounded parameter in the model allows for a better fitting for experimental curves. Indeed, we notice that, for compression times from 3 min to 5 min, the computed height/width ratio obtained for c ∼ 2 is closer to the experimental results extrapolated from [32] than the one obtained for the spherical permeability tensor (see Fig. 11a). This behaviour might be due to a better resolution of the impact of remodelling on the permeability and to the fact that, according to Eq. (45), the aggregate is "more permeable", thereby facilitating the backflow of the fluid into the aggregate. We remark, however, that the time window considered in Fig. 11a is shorter than the one in Fig. 10a because the computational effort of the simulations does not lead to an appreciable improvement of the simulated results. On the other hand, for shorter compression times (i.e. less than 3 min), for which the effects of remodelling are sufficiently weak, the introduction of the permeability tensor (45) seems to induce a decrease of the values of the pressure inside the aggregate (see Fig. 11b). Moreover, coherently with what we have said in Sect. 7.3, this might also influence the behaviour of the aggregate during the shape recovery test, thereby leading to a less spherical shape of the deformed (and remodelled) aggregate before 1 min and to a more spherical shape after 1 min.
The introduction of the permeability tensor defined in Eq. (45) affects also the pressure distribution in the aggregate, as is visible in Fig. 11b, c. In Fig. 11b, indeed, it is possible to notice that the decrease of pressure associated with the "syringe effect" is more marked than in the case of a purely spherical permeability tensor, and pressure becomes negative at an earlier time. We interpret also this result as a consequence of the fact that the aggregate is "more permeable" when Eq. (45) is employed. To be more specific, the value of the pressure at the point X U = (40,0,190) μm in the reference configuration of the aggregate becomes zero at the instant of time t = 28 s, and is negative for t > 28 s, i.e. earlier than t = 40 s (see Fig. 4).
A last remark concerns the formation of the "pressure bubble" (see Fig. 11c), which apparently shows up earlier than in the case of the spherical permeability tensor, as can be deduced by comparing Fig. 11c with Fig. 3c, both being taken at t = 30 s. As a consequence of this behaviour, the pressure values within the aggregate at t = 30 s are approximately 1/4 of those obtained with the purely spherical permeability tensor.

Concluding remarks
We addressed the mechanical response of a multicellular aggregate, modelled as a saturated and hydrated medium, featuring a porous solid phase and an interstitial fluid, and subjected to a compression-release test in displacement control. For our purposes, we considered Darcian regime for the dynamics of the fluid phase and large deformations and remodelling for the solid phase. The remodelling was identified with the anelastic phenomenon determining the evolution of the internal structure of the aggregate and culminating in the stress-driven change of its elastic parameters and in the modulation of its hydraulic properties. The analogy with plasticity offers also the possibility of investigating several theoretical branches of remodelling, which involve aspects of the problem related, for example, to Differential Geometry [24,44,102], variational methods [28,30,45,47,88] and Mathematical Analysis [29,75].
With respect to other works in the field [40,42,46], in the present contribution we proposed a threedimensional model able to solve the deformation, the interactions between the fluid and the solid phase and the contact interactions between the aggregate and the compressive apparatus. The numerical simulations were developed with the aid of the commercial software COMSOL Multiphysics TM v. 5.3a.
In response to the coupling between the plastic-like distortions associated with remodelling and the flow of the interstitial fluid, we notice the occurrence of two different, yet coupled, mechanisms. The first one appears through a "travelling pressure bubble", which characterises the circumstance in which the remodelling produces the strongest impact on the aggregate's evolution and is accompanied by a re-distribution of the stress in the aggregate itself. In fact, the "pressure bubble" distinguishes between the case in which there is a drastic reduction of the stress and the case in which, based on the results of the simulations, the mechanical properties of the aggregate are indistinguishable from those of the purely elastic case, with no structural changes. The second mechanism is the "syringe effect", which seems to descend from the interaction among the size of the aggregate, remodelling, contact and evolution of the pore space.
Another novelty of our work with respect to [42], which considered a monophasic aggregate, is that, through the visualisation of the behaviour of the interstitial fluid, we may figure out its influence on the cellular phase of the aggregate. For example, in the areas characterised by negative fluid pressure, i.e. in the proximity of the contact surfaces, one expects that the volume of the cells tends to shrink due to compression. In fact, the volume of the cells depends on the influx and/or efflux of the intracellular water and, although these aspects have not been explicitly accounted for in our model, the results of our simulations suggest to consider them as a further development of our work. Note also that the influx/efflux of intracellular water occurs over the time scale of minutes [49], which is the time scale considered in our simulations.
Looking at Figs. 7, 8 and 9, we may hypothesise the existence of interactions between the scale of the pores and that of the aggregate. To us, these are suggested by the fact that, in those figures, the mechanical response at the scale of the aggregate, expressed in terms of the force exchanged with the upper plate, is investigated for different parameters, many of which can be associated with the aggregate's microstructure. The above-mentioned interactions, in turn, influence both the final shape of the aggregate and its tendency to recover its original shape, which depends on the magnitude of stress and the re-opening of the pores (see Fig. 10).
For what concerns the stress distribution within the aggregate and the contact force exchanged between the aggregate and the plates, our results agree with those reported in [42], which employs similar geometries and material parameters for the solid phase. On the basis of this agreement, we focused also on other aspects of our model, such as on those related to the fluid, that could not be present in the monophasic setting of [42], and whose predictions could suggest new experimental investigations.
Beyond extending the poro-plastic setting, our work might suggest generalisations concerning, for example, tissue's anisotropy, other anelastic phenomena, such as damage (see e.g. [50]) or growth, and the dependence of the permeability on stress. In addition, more complex contact conditions may give rise to geometries that require gradient theories, both for the deformation and the anelastic processes, and a description of the flow capable of better resolving the behaviour of the fluid close to the boundary (Brinkman effect [14]).
As a hint for further research both for multicellular aggregates and for other hydrated biological media, we have studied the case in which the permeability tensor is an isotropic, yet non-spherical, tensor-valued function of the deformation and of the distortions associated with the remodelling. In our opinion, our results could be promising because they show how, by paying the relatively low price of a slight complication of the mathematical model, one gains two additional material parameters, i.e. k (1) 0a and k (2) 0a , that can be tuned to improve the flexibility of the model. Note that this is true even though, for our simulations, we have used only the parameter c to express k (1) 0a and k (2) 0a as a rescaling of k 0a . We conclude by saying that, although our study was successful to better approximate the experimental shape-recovery curves in the time window [0, 5] minutes, it has to be regarded as preliminary at this stage. Indeed, it is for instance still necessary to understand how the new parameters interact with the other ones. We believe, however, that our results may suggest new studies, also experimental ones, devoted to the determination of the parameters k (1) 0a and k (2) 0a that, in our present work, have been found by means of a best fit procedure.
Finally, we think that our work, especially if referred to the case of materials with a non-trivial initial distribution of stress, may be used to test the constitutive setting of [70] in which residual stresses are treated as arguments of the strain energy function (see also [83]).