Shape regulation generates elastic interaction between living cells

The organization of live cells to tissues is associated with the mechanical interaction between cells, which is mediated through their elastic environment. We model cells as spherical active force dipoles surrounded by an infinite elastic matrix, and analytically evaluate the interaction energy for different scenarios of their regulatory behavior. We obtain attraction for homeostatic (set point) forces and repulsion for homeostatic displacements. When the translational motion of the cells is regulated, the interaction energy decays with distance as $1/d^4$, while when it is not regulated the energy decays as $1/d^6$. This arises from the same reasons as the van der Waals interaction between induced electric dipoles.


A. Mechanobiological background
Live cells exert contractile forces on their environment [1]. The elastic behavior of the extracellular matrix (ECM) and its rigidity affect the forces that cells apply and consequently the resultant displacements [2][3][4][5], the cell migratory behavior [6][7][8][9][10], and cell division [11]. Alterations of the shape and size [12] and also of the rigidity [13] of cells in response to changes of the ECM rigidity were also observed experimentally. For stem cells the rigidity of the ECM may even affect their biological phenotype, from neuronal on soft substrates, muscle on substrates with intermediate rigidity to bone on rigid substrates [14]. It was also shown experimentally that the interaction behavior between cells is related to the elastic moduli and the non-linear elastic behavior of the substrate [6], and that presence of rigid boundaries affects the behavior of cells [15].
One of the immediate conclusions from these experimental results is that cells sense changes in their mechanical environment and respond to those changes in a variety of ways: by modifying the applied forces and displacements, their size and shape and even their rigidity. This, in turn suggests that there are mechanical interactions between cells through the environment [16]. Indeed, the dependence of the interaction distance between live cells on the rigidity of the substrate was observed [6]. It was also shown that on nonlinear substrates, cells respond to the presence of each other over relatively long distances (up to 10 cell diameters) [2] in contrast to linear substrates, on which the response distance is limited to 1-2 cell diameters. This fact is attributed to the strain-stiffening behavior of the ECM [2,17] or to fiber buckling [9]. Improved models of biopolymer gels predict an increase in the range of transmission of internal cellular forces [18,19]. An alternative hypothesis is that the fibrous nature of the ECM makes the transmission of mechanical signals to such long distances possible [20][21][22][23][24].
In this paper we model cells as spherical force dipoles, which create contractile forces on their surface. We seek the interaction energy between two such cells surrounded by an infinite ECM. Realistically, this ECM has nonlinear material properties, thus theoretically analyzing such interactions is a very ambitious goal [25,26]. Here we lay the conceptual foundations for reaching this goal by focusing on active force dipoles surrounded by a linearly-elastic, or Hookean material. We suggest that the concepts we introduce and the physical mechanisms that we identify would be relevant also to nonlinear media. We define the interaction energy as the additional amount of work performed by the force dipoles as a result of the presence of neighboring dipoles. From studies of the micromechanics of elastic inclusions we know that in the case of two bodies of any shape that apply a hydrostatic pressure on a surrounding linearly elastic material, their interaction energy vanishes [27,28]. We define this type of force dipoles as "dead".
As an example of such behavior one may think of heating a system with inclusions in a medium with a different thermal expansion coefficient. The self-displacements applied by each such inclusion do not depend on the distance to neighboring inclusions. In contrast to this, in this paper we introduce "live" behavior as self-regulation of the applied forces or displacements in order to preserve their shape in the presence of the interaction with other cells, or force dipoles. We calculate the interaction energy between two "live" force dipoles as a function of the separation distance between them.
The mechanical environment of the cell may change due to the presence of other cells, due to external forcing, or due to changes in the medium's rigidity, and it is not clear how cells respond to such changes. The working hypothesis of many studies, and which we will adopt here, is that there is some sort of mechanical homeostasis in the cell.
Namely, the cell tends to maintain certain quantities [29,30]. For instance, cells may regulate the forces that they apply, and then the displacement that they generate will vary with environmental changes, or alternatively, cells may regulate their deformation, and then the forces required to generate those displacements will vary. We show that shape regulation of the force dipoles leads to attractive interaction for force homeostasis and to repulsive interaction for displacement homeostasis. The interactions we find are analogous to the van der Waals interaction between two induced electric dipoles, thus corroborating the mechanobiological elasticity-electrostatics analogy [31][32][33].
It is important to emphasize that real cells probably do not keep their exact shape, however we suggest that generically self-regulation related to the interplay between active forces and the cell's shape could generate an interaction which is qualitatively similar to what we find. Our work deals with an abstract model of the cell and its mechanical behavior, which in real life are clearly more complicated. Yet we hope that insights obtained from our analytical solution of this ideal picture will shed light on the understanding of interactions in more realistic scenarios. In particular it would be interesting to relate our work to recent work on the relation between cell morphology and the polarization of the active forces that the cell applies in response to its mechanical environment [5,34,35].

B. Cells as spherical active force dipoles
In analogy to Eshelby's elastic inclusions as force-generating centers [36,37], the contractile activity of cells can be modeled as contractile force dipoles [31], see Fig. 1. We adopt this approach, and for simplicity introduce spherical force dipoles as spherical bodies, which apply isotropic contractile forces on their environment [25]. The contractile forces at the surface of each spherical force dipole represent the cellular forces that are generated by the contraction of the actin network by myosin motors and are transferred to the ECM trough focal adhesions. By analogy with electric dipoles in electrostatics, a force dipole is defined in mechanics as two equal and opposite point forces applied at some distance from each other, see Fig. 1(a). Our spherical force dipole consists of an infinite number of such linear force dipoles with the same center point, distributed isotropically on the surface of the sphere, see Fig. 1 C. No interaction between "dead" spherical force dipoles The interaction energy vanishes for two spherical inclusions in an infinite, linearly-elastic medium, each of which induces symmetric displacements or stresses in the principal directions; namely an initial volume change only with no distortion of their shapes [27,28], see Fig. 2. This result may be explained in the following way: the dilation of each sphere creates a pure shear field around it and a solely compressive field inside. The energy of interaction between the two spheres may be expressed in terms of the stress from one sphere times the strain from the other sphere, integrated over the interior of the spheres [28,38]. Since one field is a pure shear and the other is pure compression, their coupling does not generate an additional energy. Explicit evaluation of the vanishing interaction energy in this case [39] is given in Appendix A. We note that this result holds if the elastic moduli inside the inclusions are identical to those of their environment, while if the elastic moduli differ, the interaction energy will not vanish [40].
FIG. 2: Two spherical force dipoles each with radius R0, both applying a radial isotropic displacement u0 on their surfaces.
The coordinate system of sphere 1 is right handed (green) and the coordinate system of sphere 2 is left handed (blue) and may be written as θ2 = π − θ ′ 2 where θ ′ 2 is the commonly-used right-handed azimuthal coordinate for sphere 2.
For calculating the interaction energy in an infinite periodic array of spherical active force dipoles, we first note that this situation is equivalent to a single force dipole in its Wigner-Seitz unit cell, see Fig. 3. The symmetry of this periodic array dictates that the normal displacement on the midplane between each two neighboring force dipoles must vanish, and thus we may equivalently analyze a spherical force dipole surrounded by elastic ECM limited by a rigid polyhedral Wigner-Seitz unit cell. Since we describe the material as linearly elastic, if we do not take into account regulation in the activity of the force dipoles due to the deformations generated by their neighbors, we can employ the aforementioned result regarding two "dead" force dipoles and conclude that also for a system of many such force dipoles the interaction energy vanishes.

D. Spherical Wigner-Seitz unit cell approximation
Instead of solving the displacement field with the boundary condition of vanishing displacement on the actual boundaries of the polyhedral Wigner-Seitz unit cell, we previously approximated the unit cell as having a spherical shape [30]. Within this approximation we could exactly calculate the interaction energy by solving the sphericallysymmetric equations of mechanical equilibrium. We define the interaction energy as the difference between the work performed by the spherical force dipole to displace the medium by u 0 when it is surrounded by a rigid spherical Wigner-Seitz unit cell and the work when it is surrounded by an infinite ECM. For a force dipole of radius R 0 surrounded by ECM with bulk modulus K and shear modulus G. The self-energy of a single force dipole in an infinite medium is given by (see Section II E below): E 0 = 8πGR 0 u 2 0 . This is the basic energy scale in our problem, and all interaction energies we obtain scale with E 0 . For a single force dipole enclosed by a spherical Wigner-Seitz unit cell of radius R c , we relate R c to the distance between neighboring force dipoles, and define the dimensionless distance asR c = R c /R 0 . The interaction energy in the displacement-homeostasis scenario is given by [30]: is the Poisson ratio of the ECM. This interaction energy represents a repulsive force which decays algebraically with dipole-dipole separation. If instead we consider stress homeostasis and fix the stress σ 0 that the force dipole applies on the medium, we may write E 0 =

E. Motivation
At first glance it seems that these results [30] contradict the previously-mentioned general result [27,28] that bodies of any shape included in a continuous linear elastic material and applying on it isotropic dilational displacements on their surface do not interact. In order to resolve this seeming contradiction, in this paper we compare the case of a spherical force dipole inside a concentric spherical Wigner-Seitz unit cell (see Fig. 3) to the realistic configuration in which two spherical force dipoles are embedded at some distance d between them in an infinite ECM with linear properties (see Fig. 2). This comparison enables us to demonstrate that the reason for the contradiction is the difference between the definitions of the geometries in the two cases; In the two-sphere case the force dipoles obtain a drop-like shape due to the interaction [see Fig. 4(a)], while a spherical force dipole inside a spherical Wigner-Seitz unit cell preserves its spherical shape even when it interacts with other force dipoles. Since the displacements created by each force dipole distort its neighbor, shape preservation may be achieved in the two-force-dipole setup only by application of appropriate anisotropic displacements by each dipole on its surface [see Fig. 4(b),4(c)], and these give rise to the interaction energy that we study in this paper.

A. Theoretical framework
We consider two identical spheres of radius R 0 surrounded by an infinite material with linear elastic, or Hookean behavior defined by bulk modulus K and shear modulus G and denote the distance between their centers by d. In addition to the isotropic displacement u 0 , each sphere applies an anisotropic displacement, which is intended to cancel the anisotropic displacements on its surface that are caused by the other sphere. In order to simplify the calculations we choose a left-handed coordinate system for sphere 2; namely the newly defined angle θ 2 equals π − θ ′ 2 , see Fig. 2.

B. General displacements generated by spherical force dipoles
The displacement field around each force dipole must satisfy mechanical equilibrium, which we write in terms of the displacements field u as [41]: illustration purposes, the initial distance between the spheres was set to d = 3R0, and the self-displacement to u0 = 0.4R0.
The Poisson ratio is ν=0.05. We show two of the "live" cases that we describe below, which differ in whether the volumes and positions of the spheres are preserved during the interaction (b) or not (c).
Due to the symmetry with respect to rotation about the axis connecting the centers of the two force dipoles, there is no dependence on the azimuthal angle φ, thus we write Eq. (1) in spherical coordinates as: where the Laplacian in spherical coordinates excluding terms depending on φ is given by: Based on the general solution for the displacement field of a sphere with given cylindrically symmetric displacements on its surface [41], we write the anisotropic displacements field satisfying (2-3) outside the spherical force dipole (r > R 0 ) as a multipole expansion in terms of spherical harmonics Y n (θ) = 2n+1 4π P n (cos θ):  the Legendre polynomial of order n [42]. Here u ri and u θi are the radial and angular components of the displacement field caused by sphere i, and the infinite sums represent the anisotropic corrections that each force dipole produces in order to cancel the shape distortion caused by its neighbor. u 0 and R 0 have been inserted so that the coefficients C n and D n are dimensionless.
Defining the dimensionless displacements u ri = uri u0 , u θi = u θi u0 and position r = r R0 , we rewrite Eqs. (5) and (6) in dimensionless form as: Note that (7-8) solve Eq. (1) only when each force dipole is surrounded by an infinite homogeneous linearly-elastic medium, including in the interior of the neighboring force dipole. Biological cells clearly have a rigidity which differs from that of the ECM surrounding them, and thus this assumption seems to be problematic. We overcome this by realizing that we may first solve the mechanical problem in which the cells are assumed to have the same linear elastic properties as the ECM. The resultant solution includes a certain stress and displacement on the surface of each cell, and the solution outside the cells is independent on how the cell generates this stress on its surface. In particular, the stress that actual cells apply on their surrounding includes a passive stress coming from the rigidity of the cell plus an active stress coming from the external forces generated by molecular motors inside the cell. In our analysis we consider only the total stress and the work performed by it, which determines the interaction energy, and our results are valid irrespective of the mechanical rigidity of the cells themselves.

C. Cancellation condition
For our "live" force dipoles the sum of the anisotropic displacements caused by the neighbor force dipole and of all the corrections applied by the discussed force dipole must vanish on its surface. From this we derive conditions for the coefficients C n and D n so that each force dipole will preserve its spherical shape even in the presence of the interaction with its neighboring force dipole. In order to apply the cancellation condition and to derive from it the expressions for C n and D n , we transform the expressions for the displacement field of force dipole 2 to the coordinate system of force dipole 1 by substitution of the expressions for r 2 and θ 2 in terms of r 1 and θ 1 and then multiplying the displacement vector − → u 2 = (u r2 , u θ2 ) by the rotation matrix: We then write the resultant expressions for the radial and angular displacements caused by force dipole 2 on the surface of force dipole 1 in terms of the spherical harmonics of sphere 1 by writing: As may be seen from (10,11), every spherical-harmonic mode of force dipole 2 contributes to all the modes on the surface of force dipole 1.
As explained above, the sum of the anisotropic displacements caused by force dipole 2 must be canceled on the surface of force dipole 1 by the corrections that it applies. We write the dimensionless displacement u 11 created by force dipole 1 on its surface (namely at r 1 = 1): The term δ n,0 in (12) is a Kronecker delta, which represents the isotropic radial displacement created by sphere 1 on its surface without the anisotropic cancellation corrections. This constant term does not depend on changes in the environment of the force dipole. The remaining terms are different modes of additional displacement created by this "live" force dipole in response to the displacement field induced on its surface by the neighboring force dipole. The dimensionless displacement u 21 created by force dipole 2 on the surface of force dipole 1 is: where the sum over m originates from the fact that the displacement u 2 created by force dipole 2 is given by a multipole expansion (7)(8) with the corrective magnitudes C m and D m . The sum over n originates from the fact that after the coordinate transformation, when these modes are expressed in terms of the spherical harmonics in the coordinate system of force dipole 1, each mode from force dipole 2 contributes to all the modes of force dipole 1.
between the spherical force dipoles.
We now require that for "live" force dipoles the total displacement u 11 (θ 1 ) + u 21 (θ 1 ) on the surface of force dipole 1 is isotropic. We begin by considering the simplest (but strictest) regulation scenario, for which not only is this total displacement isotropic, but its magnitude remains equal to the displacement u 0 in the absence of interactions between the force dipoles. Moreover, we require that the center of symmetry of each force dipole does not move. We will later consider three additional scenarios in which the interaction causes the force dipoles to change their volume and/or to move, yet they remain spherically symmetric. Thus at this point we require that: Substituting (12,13,14,15) in (16,17) yields: Due to the orthogonality of the Legendre polynomials, for these infinite sums to satisfy these cancellation conditions, each term in the sums must cancel independently. Thus for all n ≥ 1 we require: Note that for n = 0, from (7-8) C 0 is irrelevant, thus we set it to zero. Moreover, since Y 0 (θ 1 ) = 1, dY0(θ1) dθ1 = 0 and (19) holds trivially, thus for n = 0 we obtain only one equation, from Eq. (18): We obtain closure of the infinite coupled linear Eqs. (20,21,22) by assuming that C n = 0 and D n = 0 for n > n max .
This is justified since we will be interested in large separations between the force dipoles, and due to the fact that the solutions decay as 1/r n , thus at large r, large n terms become negligible. We verify this numerically by increasing n max until convergence, see Section III below. to have an additional contribution due to the interaction. Similarly we consider either fixed position (FP) or variable position (VP) of the force dipole (n = 1). In the VP case instead of the cancellation condition (20)(21), the coefficients C 1 , D 1 obey the following: Equations (23-24) allow us to find C 1 and D 1 that cancel the deformation at mode 1 without canceling the translational motion of the force dipole. They follow from Eqs. (20)(21) and from the fact that for translation the coefficients of P 1 (cosθ) = cos(θ) in the radial direction and of dP1(cosθ) dθ = sin(θ) in the tangential direction must be equal.

E. Interaction energy
After obtaining the coefficients C n and D n and substituting them in Eqs. (12)(13)(14)(15) we compute the interaction energy from the work performed by the force dipoles to generate this deformation: Here the integration is over the surfaces of both spheres, − → F is the force per unit area applied by each force dipole on its environment. It is important to emphasize the difference in calculation of the interaction energies in the cases of dead versus live force dipoles. In the general case of two dead force dipoles the displacements and forces applied by each one of them are not affected by changes in its environment and thus the amount of additional work done by it equals [43] [see Eq. (25)]: In the particular case that the force dipoles apply isotropic forces or displacements causing only volume change without distortion (as for two "dead" spherical force dipoles) the interaction energy ∆E vanishes [27]. In contrast, for two "live" force dipoles the forces and displacements created by each one of them are modified due to the interaction with the neighbor. The resultant forces and displacements applied by "live" force dipoles are not isotropic and the interaction energy does not vanish anymore. Here the additional work that is performed by each force dipole is: Here δ − → u and δ − → F are the "live" parts of the displacement and the force created by the force dipoles due to their interaction. For two "live" force dipoles, we evaluate the interaction energy using the stress tensor τ that arises in In the FS cases D 0 is evaluated from the cancellation condition and Eq. (28) becomes: while in the VS cases D 0 is zero and Eq. (28) becomes: The meaning of the sum of the last two terms in Eqs. (29)(30) is the change of the volume of force dipole 1 caused by all the modes of displacement created by force dipole 2 (see Fig. 2). We thus define: Thus the interaction energies in the FS cases become: and in the VS cases:

III. RESULTS
We evaluated the interaction energy by terminating the infinite sums (20)(21) at different values of n max , up to n max = 15 for all four regulation scenarios. Figure 6 shows the normalized interaction energy ∆Ẽ = ∆E/E 0 vs the normalized distance between the force dipolesd = d/R 0 for all four regulation scenarios. Figure 7 shows the convergence of the interaction energy with increasing the value of n max . Even at the smallest dipole-dipole separations the solution converges with moderate value of n max = 4 − 5. Moreover, for distances larger than 4 sphere radii in the FSVP, VSFP and VSVP cases and for distances larger than 6 sphere radii in the FSFP case the energy is well approximated by taking the minimal n max possible. Namely for the FP cases n max = 1 suffices, and for the VP cases we take n max = 2, see Table I Excluding the n = 1 term of the sum in Eqs. (5)(6) constitutes the VP cases in which the force dipoles are free to move but preserve their size and shape, and excluding the n = 0 term constitutes the VS cases in which the force dipoles preserve their spherical shape but are free to change their volume due to the interaction. We see from Eq. (28) that the presence of the coefficient D 0 determines the sign of the interaction energy, see also Fig. 6 and Table I This is similar to our analysis within the spherical unit-cell approximation of different homeostasis scenarios, where we found repulsion when the force dipole generates a fixed displacement on its boundary and attraction when it applies a fixed stress on its environment [30]. Relating this to our present results, in the FS cases the resultant displacement on the boundary of each force dipole is fixed and does not depend on the changes in its environment, which corresponds to displacement homeostasis. In the VS cases the size of each force dipole changes due to the interaction, namely the displacement is not fixed. There is no correction to the n = 0 mode, which means that there is no additional active force or stress at this mode, and we relate this to stress homeostasis, in which indeed the interaction is attractive.
is the dimensionless distance between the force dipoles, and ∆V21 is defined in Eq. (33).
The corresponding expressions for arbitrary distance are given in Appendix E.
An intuitive explanation for the sign difference in the interaction may be considered as follows; FS means that in the proximity to other force dipoles, each force dipole has to exert an additional force in order to generate the displacement that it was programmed to have. Thus more energy is required and the force dipoles would benefit energetically from moving away from each other. As noted above, VS is related to stress homeostasis, for which the tension that neighboring contractile force dipoles generate around themselves add up, and less energy is required in order to reach the homeostatic stress value, thus it is beneficial for cells to be close to other cells.
The position regulation (n = 1) sets the strength of these attractive or repulsive interactions. For VP, C 1 = D 1 = 0 and at long distances the interaction energy decays as 1/d 6 , while for FP, C 1 and D 1 are nonzero, which leads to a 1/d 4 decay, see Fig. 6 and Table I. To quantitatively relate these functional forms to the results of our spherical unit-cell approximation, we consider an infinite array of active force dipoles. Using our results for the interaction between two active force dipoles, we now evaluate the interaction energy in an infinite 3D array of such shape-regulating spherical force dipoles (see Fig. 3). In a linear medium and for small displacements generated by each force dipole, by superposition the total interaction energy per force dipole is: where n is the density of force dipoles at distance r from any given force dipole, which for large distances we assume is uniform, and we introduce a near-field cutoff distance equal to the lattice spacing D.
Substituting the expressions for the two-force-dipole interactions from Table I  The expressions for the interaction energy given in Table I (and similarly in Table II in Appendix E for arbitrary distance) all vanish in the incompressible limit (ν = 0.5). The self-displacement generated by each force dipole separately includes a pure compressive mode inside that force dipole, and sustaining this in an incompressible medium with a fixed displacement u 0 on the boundary requires an infinite force. Thus the problem of mechanical interactions between force dipoles in a completely incompressible medium deserves a separate analysis, which we defer to future work.

IV. DISCUSSION
We model live cells as spherical force dipoles surrounded by a linear, or Hookean elastic environment. Since in the case of isotropic displacements and forces the interaction energy between force dipoles vanishes, we distinguish between "dead" behavior in which the force dipoles apply constant forces and self-displacements on their surface, and "live" behavior in which the forces and self-displacements applied change in response to changes in the environment of the force dipoles. We solved the interaction energy for four different types of such regulation in which the force dipoles preserve their spherical shape and in addition volume, position, both or none of them. We found the interaction energy to be inversely proportional to the distance between the force dipoles to the fourth power in the case of fixed position and to the sixth power in the case of variable position. These results are similar to and stem from the same reasons as the van der Waals dipole-induced dipole interaction in electrostatics [44]. We also found that in the case of volume preservation the force dipoles are repelled while without volume preservation they are attracted.
We solved the deformation fields for the case when the rigidity of the force dipoles is identical to that of their environment. However, biological cells are complex entities and their rigidity differs from point to point and also differs from the rigidity of the extracellular matrix. In order to relate our results to live cells we describe each of them as a mechanism that applies forces on its environment on the surface and responds by their variation to application of external force or displacement. The displacements and forces applied by a cell may be divided into "passive" and "active" parts. The "passive" part of the forces or displacements would stay the same if the cells were dead and preserve their elastic properties, while the "active" part depends on the programmed behavior of the cells and is generated by the contraction of acto-myosin networks inside them. Since the resultant force and displacement are the sum of those two parts cells may create such "active" response such that the resultant forces and displacements will coincide with the case considered here, for which their rigidity coincides with the rigidity of their environment.
Stress homeostasis vs displacement homeostasis is believed to occur in different cell types and in different mechanical environments [45]. We suggest to examine in experiments the correlation that we find here between these distinct regulatory behaviors and attraction vs repulsion. Specifically, for cells that are known to regulate their stress or displacement one could test whether they indeed are attracted or repelled, respectively. Alternatively, when the regulatory behavior is not known, we suggest that our results will enable to infer it from the direction of the interaction.
It would be interesting to consider also more complicated regulatory behaviors on top of the two extreme limits of fixed stress and fixed displacement. In particular, it would be interesting to test whether our fixed-position vs variableposition cases could generate more complex scenarios, since we predict that that would relate to the interaction strength. Clearly the interaction energy that we focus our analysis on may not be directly probed in experiments. We suggest that our predictions on attraction or repulsion between cells could be tested by studying other experimental indications of inter-cellular interactions. For instance, cells that are attracted to each other tend to try and move closer together but in rigid environments may not move and instead send protrusions one toward each other [6].
Moreover, focusing of mechanical stress between interacting cells may be experimentally visualized [2].
Following our work on three-dimensional spherical force dipoles, it would be interesting to solve the case of twodimensional disks on an elastic medium. The same four regulatory possibilities can be considered. This could more directly be related to experiments in which cells are grown on surfaces of different elastic materials [6]. We expect that in this case the interaction energy will not vanish even for "dead" force dipoles. Clearly the work presented in this paper is a first step toward theoretical analysis of interactions that stem from mechanobiological regulation, and in particular shape regulation of live cells. Our focus on linearly-elastic, homogeneous and isotropic media enabled us to obtain the detailed analytical description presented above. It would be important to continue this line of research and study the influence of the material nonlinearity, heterogeneity and anisotropy of the ECM on the interaction between cells.
In order to evaluate the resultant displacement of two such force dipoles we superimpose the displacement fields of both force dipoles. To do so, we translate the displacement field created by force dipole 1 to the coordinate system of force dipole 2. Using the cosine theorem we rewrite the displacement fields created by the force dipoles in terms of θ 2 and d, see Fig. 2: Here r 1 and r 2 are unit vectors in the radial directions of force dipoles 1 and 2, respectively. The resultant displacement in the radial direction of force dipole 2 is:

20
The symmetry of the system imposes that the total work done by the system reads: After evaluation we obtain: Thus the interaction energy vanishes in this case. A similar analysis may be done for two force dipoles which apply radial displacements u 0 on their surfaces. The energy in this case may be written as: and also here the interaction energy vanishes.
Using (23)(24) we get in these cases: The last two terms in Eq. (D4) are proportional to higher powers of C n and D n . It may be seen from Table I that the coefficients C n and D n are inversely proportional to d, and thus at large distances the last two terms in Eq. (D4) may be neglected and Eq. (D4) reduces to Eq. (D3).
Appendix E: D0, C1, D1, C2, D2 and ∆E for arbitrary distance In Table II we present the expressions for the coefficients D 0 , C 1 , D 1 , C 2 , D 2 and for the dimensionless interaction energy ∆ E without the large-distance assumption that appears in Table I.  into account, ∆ E ≡ E int E 0 is the resultant dimensionless interaction energy and ∆ E∞ is its asymptotic long-distance behavior given also in Table I. d = d R 0 is the dimensionless distance between the force dipoles.
We use the following expressions: