Shape regulation generates elastic interaction between living cells

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–5], the cell migratory behavior [6–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 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 has also been 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. Such sensing, response and regulation is performed at the microscopic (sub-micrometer) level by the actomyosin cytoskeleton, which generates variable contractile stresses [16–18], and furthermore by the adhesion of cells to their environment [19, 20]. Macro-scale (tens of micrometers) experiments suggest that there are mechanical interactions between cells through the environment [21]. Indeed, dependence of the interaction distance between live cells on the rigidity of the substrate has been observed [6]. It has also been 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, 22] or to fiber buckling [9]. Improved models of biopolymer gels predict an increase in the range of transmission of internal cellular forces [23, 24]. An alternative hypothesis is that the fibrous nature of the ECM makes the transmission of mechanical signals over such long distances possible [25–29].

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 [30, 31]. 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 should also be relevant 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 [32, 33]. We define this type of force dipole 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 and of their regulatory behavior, as detailed below.

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. That is to say, the cell tends to maintain certain quantities [34, 35]. 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 [36–38].

It is important to emphasize that real cells may modulate both their shape and their contractility in response to changes in their environment such as changes in the rigidity of the extracellular matrix or application of external forces on the cells. Constructing models that describe such morphological changes requires analyzing the internal structure of the cells. Instead, 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 cells apply in response to their mechanical environment [5, 39–42].

In analogy to Eshelby's elastic inclusions as force-generating centers [43, 44], the contractile activity of cells can be modeled as contractile force dipoles [36], see figure 1. We adopt this approach, and for simplicity introduce spherical force dipoles as spherical bodies, which apply isotropic contractile forces on their environment [30]. 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 figure 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 figure 1(b).

Zoom In Zoom Out Reset image size

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; i.e. an initial volume change only with no distortion of their shapes [32, 33], see figure 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 all space [33, 45]. Since one field is a pure shear and the other is pure compression, their coupling does not generate any additional energy. Explicit evaluation of the vanishing interaction energy in this case [46] 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 [47].

Zoom In Zoom Out Reset image size

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 figure 3. The symmetry of this periodic array dictates that the normal displacement on the midplane between each two neighboring force dipoles must vanish. Thus we may equivalently analyze a spherical force dipole surrounded by elastic ECM, which in turn is 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 the interaction energy also vanishes for a system of many such force dipoles.

Zoom In Zoom Out Reset image size

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 [35]. Within this approximation we could exactly calculate the interaction energy by solving the spherically-symmetric 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₀ 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₀ surrounded by ECM with shear modulus G, the self-energy of a single force dipole in an infinite medium is given by (see section 2.5 below): ${E}_{0}=8\pi {{GR}}_{0}{u}_{0}^{2}$. This is the basic energy scale in our problem, and all interaction energies we obtain scale with E₀. 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 as ${\tilde{R}}_{{\rm{c}}}={R}_{{\rm{c}}}/{R}_{0}$. The interaction energy in the displacement-homeostasis scenario is given by [35]: ${\rm{\Delta }}E=\tfrac{3-\nu }{4(1-\nu )}\tfrac{{E}_{0}}{{\widetilde{R}}_{{\rm{c}}}^{3}}$, where $\nu =\tfrac{3K-2G}{2(3K+G)}$ is the Poisson ratio and K the bulk modulus 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 ${\sigma }_{0}$ that the force dipole applies on the medium, we may write ${E}_{0}=\tfrac{\pi {\sigma }_{0}^{2}{R}_{0}^{3}}{2G}$ and now ${\rm{\Delta }}E=-\tfrac{3-\nu }{4(1-\nu )}\tfrac{{E}_{0}}{{\widetilde{R}}_{{\rm{c}}}^{3}}$, representing an attractive interaction.

At first glance it seems that these results [35] contradict the previously-mentioned general result [32, 33] 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 figure 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 figure 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 figure 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 figures 4(b), (c)), and these give rise to the interaction energy that we study in this paper.

Zoom In Zoom Out Reset image size

We consider two identical spheres of radius R₀ surrounded by an infinite material with linear elastic (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₀, 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: specifically, the newly defined angle ${\theta }_{2}$ equals $\pi -{\theta }_{2}^{\prime }$, see figure 2.

The displacement field around each force dipole must satisfy mechanical equilibrium, which we write in terms of the displacement field $\vec{u}$ as [48]:

$$\begin{eqnarray}&&\displaystyle \frac{1}{1-2\nu }{\rm{\nabla }}{\rm{\nabla }}\cdot \overrightarrow{u}+{{\rm{\nabla }}}^{2}\overrightarrow{u}=0.\end{eqnarray} \tag{ 1 }$$

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 equation (1) in spherical coordinates as:

$$\begin{eqnarray}&&\displaystyle \frac{1}{1-2\nu }\displaystyle \frac{\partial }{\partial r}\left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}({r}^{2}{u}_{r})+\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial }{\partial \theta }({u}_{\theta }\sin \theta )\right]+{{\rm{\nabla }}}^{2}{u}_{r}-\displaystyle \frac{2}{{r}^{2}}{u}_{r}-\displaystyle \frac{2}{{r}^{2}}\displaystyle \frac{\partial {u}_{\theta }}{\partial \theta }-\displaystyle \frac{2{u}_{\theta }\cot \theta }{{r}^{2}}=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{1-2\nu }\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial \theta }\left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}({r}^{2}{u}_{r})+\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial }{\partial \theta }({u}_{\theta }\sin \theta )\right]+{{\rm{\nabla }}}^{2}{u}_{\theta }+\displaystyle \frac{2}{{r}^{2}}\displaystyle \frac{\partial {u}_{r}}{\partial \theta }-\displaystyle \frac{{u}_{\theta }}{{r}^{2}{\sin }^{2}\theta }=0,\end{eqnarray} \tag{ 3 }$$

where the Laplacian in spherical coordinates excluding terms depending on ϕ is given by:

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}=\displaystyle \frac{1}{{r}^{2}\sin \theta }\left[\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\sin \theta \displaystyle \frac{\partial }{\partial r}\right)+\displaystyle \frac{\partial }{\partial \theta }\left(\sin \theta \displaystyle \frac{\partial }{\partial \theta }\right)\right].\end{eqnarray} \tag{ 4 }$$

Based on the general solution for the displacement field of a sphere with given cylindrically-symmetric displacements on its surface [48], we write the anisotropic displacements field satisfying (2)–(3) outside the spherical force dipole ($r\gt {R}_{0}$) as a multipole expansion in terms of spherical harmonics ${Y}_{n}(\theta )=\sqrt{\tfrac{2n+1}{4\pi }}{P}_{n}(\cos \theta )$:

$$\begin{eqnarray}&&{u}_{{ri}}=\displaystyle \frac{{u}_{0}{R}_{0}^{2}}{{r}_{i}^{2}}+{u}_{0}\displaystyle \sum _{n=0}^{\infty }\left[n(n+3-4\nu )\displaystyle \frac{{C}_{n}{R}_{0}^{n}}{{r}_{i}^{n}}-(n+1)\displaystyle \frac{{D}_{n}{R}_{0}^{n+2}}{{r}_{i}^{n+2}}\right]{Y}_{n}({\theta }_{i}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{u}_{\theta i}={u}_{0}\displaystyle \sum _{n=0}^{\infty }\left[(-n+4-4\nu )\displaystyle \frac{{C}_{n}{R}_{0}^{n}}{{r}_{i}^{n}}+\displaystyle \frac{{D}_{n}{R}_{0}^{n+2}}{{r}_{i}^{n+2}}\right]\displaystyle \frac{{{\rm{d}}{Y}}_{n}({\theta }_{i})}{{\rm{d}}{\theta }_{i}},\end{eqnarray} \tag{ 6 }$$

with

$$\begin{eqnarray}&&{P}_{n}(x)={2}^{n}\cdot \displaystyle \sum _{k=0}^{n}{x}^{k}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}\tfrac{n+k-1}{2}\\ n\end{array}\right)\end{eqnarray} \tag{ 7 }$$

the Legendre polynomial of order n [49]. Here ${u}_{{ri}}$ and ${u}_{\theta i}$ are the radial and angular components of the displacement field caused by force dipole 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₀ and R₀ have been inserted so that the coefficients C_n and D_n are dimensionless.

Defining the dimensionless displacements ${\widetilde{u}}_{{ri}}=\tfrac{{u}_{{ri}}}{{u}_{0}},{\widetilde{u}}_{\theta i}=\tfrac{{u}_{\theta i}}{{u}_{0}}$ and position $\widetilde{r}=\tfrac{r}{{R}_{0}}$, we rewrite equations (5) and (6) in dimensionless form as:

$$\begin{eqnarray}&&{\widetilde{u}}_{{ri}}({\tilde{r}}_{i},{\theta }_{i})=\displaystyle \frac{1}{{\widetilde{r}}_{i}^{2}}+\displaystyle \sum _{n=0}^{\infty }\left[n(n+3-4\nu )\displaystyle \frac{{C}_{n}}{{\widetilde{r}}_{i}^{n}}-(n+1)\displaystyle \frac{{D}_{n}}{{\tilde{r}}_{i}^{n+2}}\right]{Y}_{n}({\theta }_{i}),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\widetilde{u}}_{\theta i}({\tilde{r}}_{i},{\theta }_{i})=\displaystyle \sum _{n=0}^{\infty }\left[(-n+4-4\nu )\displaystyle \frac{{C}_{n}}{{\widetilde{r}}_{i}^{n}}+\displaystyle \frac{{D}_{n}}{{\tilde{r}}_{i}^{n+2}}\right]\displaystyle \frac{{{dY}}_{n}({\theta }_{i})}{d\theta }.\end{eqnarray} \tag{ 9 }$$

Note that (8)–(9) solve equation (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 surroundings 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.

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₂ and ${\theta }_{2}$ in terms of r₁ and ${\theta }_{1}$ followed by multiplication of the displacement vector $\overrightarrow{{u}_{2}}=({u}_{r2},{u}_{\theta 2})$ by the rotation matrix:

$$\begin{eqnarray}{\bf{B}}=\left(\begin{array}{cc}-\,\cos ({\theta }_{1}+{\theta }_{2}) & \sin ({\theta }_{1}+{\theta }_{2})\\ \sin ({\theta }_{1}+{\theta }_{2}) & \cos ({\theta }_{1}+{\theta }_{2})\end{array}\right).\end{eqnarray} \tag{ 10 }$$

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 as:

$$\begin{eqnarray}&&{({u}_{r})}_{n}=2\pi {\int }_{0}^{\pi }{u}_{r}(\theta ){Y}_{n}(\theta )\sin \theta {\rm{d}}\theta ,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{({u}_{\theta })}_{n}=\displaystyle \frac{2\pi }{n(n+1)}{\int }_{0}^{\pi }{u}_{\theta }(\theta )\displaystyle \frac{{Y}_{n}(\theta )}{{\rm{d}}\theta }\sin \theta {\rm{d}}\theta .\end{eqnarray} \tag{ 12 }$$

As may be seen from (11), (12), 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 ${\widetilde{u}}_{11}$ created by force dipole 1 on its surface (specifically at ${\widetilde{r}}_{1}=1$):

$$\begin{eqnarray}&&{\widetilde{u}}_{r11}({\theta }_{1})\equiv {\widetilde{u}}_{r1}(1,{\theta }_{1})=\displaystyle \sum _{n=0}^{\infty }[n(n+3-4\nu ){C}_{n}-(n+1)({D}_{n}-\sqrt{4\pi }{\delta }_{n,0})]{Y}_{n}({\theta }_{1}),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\widetilde{u}}_{\theta 11}({\theta }_{1})\equiv {\widetilde{u}}_{\theta 1}(1,{\theta }_{1})=\displaystyle \sum _{n=0}^{\infty }[(-n+4-4\nu ){C}_{n}+{D}_{n}]\displaystyle \frac{{{\rm{dY}}}_{n}({\theta }_{i})}{{\rm{d}}\theta }.\end{eqnarray} \tag{ 14 }$$

The term ${\delta }_{n,0}$ in (13) 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 ${\widetilde{u}}_{21}$ created by force dipole 2 on the surface of force dipole 1 is:

$$\begin{eqnarray}&&{\widetilde{u}}_{r21}({\theta }_{1})=\displaystyle \sum _{n=0}^{\infty }\,\displaystyle \sum _{m=0}^{\infty }[{f}_{{nm}}^{{Cr}}(\widetilde{d}){C}_{m}+{f}_{{nm}}^{{Dr}}(\widetilde{d})({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})]{Y}_{n}({\theta }_{1}),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\widetilde{u}}_{\theta 21}({\theta }_{1})=\displaystyle \sum _{n=0}^{{\rm{\infty }}}\,\displaystyle \sum _{m=0}^{{\rm{\infty }}}[{f}_{nm}^{C\theta }(\widetilde{d}){C}_{m}+{f}_{nm}^{D\theta }(\widetilde{d}){D}_{m}]\displaystyle \frac{{{\rm{d}}Y}_{n}({\theta }_{1})}{{\rm{d}}{\theta }_{1}},\end{eqnarray} \tag{ 16 }$$

where the sum over m originates from the fact that the displacement ${\widetilde{u}}_{2}$ created by force dipole 2 is given by a multipole expansion (8)–(9) 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. The functions ${f}_{{nm}}^{{Cr}}(\widetilde{d}),{f}_{{nm}}^{{Dr}}(\widetilde{d}),{f}_{{nm}}^{C\theta }(\widetilde{d})$ and ${f}_{{nm}}^{D\theta }(\widetilde{d})$, given in appendix B depend only on the dimensionless distance $\widetilde{d}=\tfrac{d}{{R}_{0}}$ between the two spherical force dipoles.

We now require that for 'live' force dipoles the total displacement ${\widetilde{u}}_{11}({\theta }_{1})+{\widetilde{u}}_{21}({\theta }_{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₀ in the absence of interactions between the force dipoles. Moreover, we require that the center of symmetry of each force dipole does not move. This will be denoted fixed-size, fixed-position. 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, see figure 5.

Zoom In Zoom Out Reset image size

Thus at this point we require that:

$$\begin{eqnarray}&&{\widetilde{u}}_{r11}({\theta }_{1})+{\widetilde{u}}_{r21}({\theta }_{1})\equiv 1,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\widetilde{u}}_{\theta 11}({\theta }_{1})+{\widetilde{u}}_{\theta 21}({\theta }_{1})\equiv 0.\end{eqnarray} \tag{ 18 }$$

Substituting (13)–(16) in (17), (18) yields:

$$\begin{eqnarray}&&\displaystyle \sum _{n=0}^{\infty }\{[n(n+3-4\nu ){C}_{n}-(n+1)({D}_{n}-\sqrt{4\pi }{\delta }_{n,0})]\\ &&\quad \left.+\displaystyle \sum _{m=0}^{\infty }[{f}_{{nm}}^{{Cr}}(\widetilde{d}){C}_{m}+{f}_{{nm}}^{{Dr}}(\widetilde{d})({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})]\right\}{Y}_{n}({\theta }_{1})=1,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{n=0}^{{\rm{\infty }}}\left\{[(-n+4-4\nu ){C}_{n}+{D}_{n}]+\displaystyle \sum _{m=0}^{{\rm{\infty }}}[{f}_{nm}^{C\theta }(\widetilde{d}){C}_{m}+{f}_{nm}^{D\theta }(\widetilde{d}){D}_{m}]\right\}\displaystyle \frac{{{\rm{d}}Y}_{n}({\theta }_{1})}{{\rm{d}}{\theta }_{1}}=0.\end{eqnarray} \tag{ 20 }$$

Due to the orthogonality of the Legendre polynomials, for these infinite sums to satisfy the cancellation conditions, each term in the sums must cancel independently. Thus for all $n\geqslant 1$ we require:

$$\begin{eqnarray}&&n(n+3-4\nu ){C}_{n}-(n+1){D}_{n}+\displaystyle \sum _{m=0}^{\infty }[{f}_{{nm}}^{{Cr}}{C}_{m}+{f}_{{nm}}^{{Dr}}({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})]=0,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&(-n+4-4\nu ){C}_{n}+{D}_{n}+\displaystyle \sum _{m=0}^{\infty }[{f}_{{nm}}^{C\theta }{C}_{m}+{f}_{{nm}}^{D\theta }{D}_{m}]=0.\end{eqnarray} \tag{ 22 }$$

Note that for n = 0, from (8)–(9) C₀ is irrelevant, thus we set it to zero. Moreover, since ${Y}_{0}({\theta }_{1})=1,\tfrac{{{\rm{d}}{Y}}_{0}({\theta }_{1})}{{\rm{d}}{\theta }_{1}}=0$ and (20) holds trivially, thus for n = 0 we obtain only one equation, from equation (19):

$$\begin{eqnarray}&&-{D}_{0}+\displaystyle \sum _{m=0}^{\infty }[{f}_{0m}^{{Cr}}{C}_{m}+{f}_{0m}^{{Dr}}({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})]=0.\end{eqnarray} \tag{ 23 }$$

We obtain closure of the infinite coupled linear equations (21)–(23) by assuming that C_n = 0 and D_n = 0 for $n\gt {n}_{\max }$, with some arbitrary value of n_max, which will determine the accuracy of our calculation. 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 3 below. We evaluate the forces created by the force dipole using equations (C1)–(C6) in appendix C for the stress tensor in the elastic environment. Due to force balance, the active force per unit area generated by the force dipole is exactly minus this elastic stress.

We emphasize the physical meaning of the first two terms in our solution. Each of the spherical harmonics ${Y}_{n}(\theta )$ represents a different mode of deformation of the force dipole. For scalar fields the first two modes represent volume change (n = 0) and translation of the center of symmetry (n = 1). However, since in our case the displacement field is vectorial, the n = 1 mode may represent a combination of translation and deformation. Our definition of 'live' force dipoles focuses only on shape regulation, thus inclusion in the solution of the terms n = 0 and the part of the n = 1 term that controls translation is not mandatory. Inclusion of one or both of these terms will result in complementary condition(s) in addition to the shape preservation condition, and will describe a different type of self-regulation of the force dipoles. Overall, there are four combinations for the inclusion or exclusion of the n = 0 and n = 1 terms, and thus we analyze four possible scenarios of homeostasis or self-regulation—see figures 4(b)–(c) and 5 above.

Exclusion of one of the terms from the solution means that the corresponding correction terms C_n and D_n vanish, and thus equations (21)–(23) will not vanish anymore at the corresponding modes, i.e. non-zero resultant displacement will exist (see below). At the same time, the force applied by the dipole and the resulting elastic stress will vanish at these modes. Mode n = 0 of the displacement never vanishes, since it includes the constant term u₀ representing the fixed symmetric displacement created by the force dipole, which does not depend on the interaction.

For the size $(n=0)$ we have two options—fixed size (FS), i.e. displacement homeostasis, in which we regulate this mode, and variable size (VS) i.e. stress homeostasis, in which we allow the displacement 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 FS case the coefficient D₀ obeys the corresponding cancellation condition (23), while in the VS case it vanishes (${C}_{0}=0$ in both cases). In the VP case, instead of the cancellation conditions (21)–(22), the coefficients C₁ and D₁ obey the following:

$$\begin{eqnarray}&&(4-4\nu ){C}_{1}-2{D}_{1}+\displaystyle \sum _{m=0}^{\infty }[{f}_{1m}^{{Cr}}{C}_{m}+{f}_{1m}^{{Dr}}({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})]=\displaystyle \sum _{m=0}^{\infty }\left(\displaystyle \frac{{f}_{1m}^{{Cr}}+{f}_{1m}^{C\theta }}{2}{C}_{m}+\displaystyle \frac{{f}_{1m}^{{Dr}}+{f}_{1m}^{D\theta }}{2}{D}_{m}\right),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&(3-4\nu ){C}_{1}+{D}_{1}+\displaystyle \sum _{m=0}^{\infty }[{f}_{1m}^{C\theta }{C}_{m}+{f}_{1m}^{D\theta }{D}_{m}]=\displaystyle \sum _{m=0}^{\infty }\left(\displaystyle \frac{{f}_{1m}^{{Cr}}+{f}_{1m}^{C\theta }}{2}{C}_{m}+\displaystyle \frac{{f}_{1m}^{{Dr}}+{f}_{1m}^{D\theta }}{2}{D}_{m}\right).\end{eqnarray} \tag{ 25 }$$

Equations (24)–(25) allow us to find C₁ and D₁ that cancel the deformation at mode 1 without canceling the translational motion of the force dipole. They follow from equations (21)–(22) and from the fact that for translation the coefficients of ${P}_{1}(\cos \theta )=\cos (\theta )$ in the radial direction and of $\tfrac{{{\rm{d}}{P}}_{1}(\cos \theta )}{{\rm{d}}\theta }=-\sin (\theta )$ in the tangential direction must be equal.

After finding the coefficients C_n and D_n, we may substitute them in equations (13), (15) to obtain the n = 0 component of the total displacement

$$\begin{eqnarray}&&g=\displaystyle \frac{1}{\sqrt{4\pi }}\displaystyle \sum _{m=0}^{\infty }[{f}_{0m}^{{Cr}}{C}_{m}+{f}_{0m}^{{Dr}}({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})],\end{eqnarray} \tag{ 26 }$$

and that of the n = 1 mode

$$\begin{eqnarray}&&h=\sqrt{\displaystyle \frac{3}{4\pi }}\left\{(4-4\nu ){C}_{1}-2{D}_{1}+\displaystyle \sum _{m=0}^{\infty }[{f}_{1m}^{{Cr}}{C}_{m}+{f}_{1m}^{{Dr}}({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})].\right\}\end{eqnarray} \tag{ 27 }$$

Namely, in the VS cases the total symmetric radial displacement is $(1+g){u}_{0}$, and in the VP cases the center of symmetry is shifted by hu₀.

After obtaining the coefficients C_n and D_n and substituting them in equations (13)–(16) we compute the interaction energy from the extra work performed by the force dipoles to generate this deformation:

$$\begin{eqnarray}&&{\rm{\Delta }}E=E-2{E}_{0}=\displaystyle \frac{1}{2}{\int }_{S}(\overrightarrow{u}\cdot \overrightarrow{F}-\overrightarrow{{u}_{0}}\cdot \overrightarrow{{F}_{0}}){ds}.\end{eqnarray} \tag{ 28 }$$

Here the integration is over the surfaces of both spheres, $\overrightarrow{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. For 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 [50] (see equation (28)):

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{i}=\displaystyle \frac{1}{2}{\int }_{{S}_{i}}[(\overrightarrow{{u}_{i}}+\overrightarrow{{u}_{j}})\cdot \overrightarrow{{F}_{i}}-\overrightarrow{{u}_{i}}\cdot \overrightarrow{{F}_{i}}]{{\rm{d}}{S}}_{i}=\displaystyle \frac{1}{2}{\int }_{{S}_{i}}\overrightarrow{{u}_{j}}\cdot \overrightarrow{{F}_{i}}{{\rm{d}}{S}}_{i}.\end{eqnarray} \tag{ 29 }$$

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 ${\rm{\Delta }}E$ vanishes [32]. 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 no longer vanishes. Here the additional work performed by each force dipole is:

$$\begin{eqnarray}{\rm{\Delta }}{E}_{i} & = & \displaystyle \frac{1}{2}{\displaystyle \int }_{{S}_{i}}[(\overrightarrow{{u}_{i}}+\delta \overrightarrow{{u}_{i}}+\overrightarrow{{u}_{j}}+\delta \overrightarrow{{u}_{j}})(\overrightarrow{{F}_{i}}+\delta \overrightarrow{{F}_{i}})-\overrightarrow{{u}_{i}}\cdot \overrightarrow{{F}_{i}}]{{\rm{d}}{S}}_{i}\\ & = & \displaystyle \frac{1}{2}{\displaystyle \int }_{{S}_{i}}[\overrightarrow{{u}_{j}}\cdot \overrightarrow{{F}_{i}}+\delta \overrightarrow{{u}_{i}}\cdot \overrightarrow{{F}_{i}}+\delta \overrightarrow{{u}_{j}}\cdot \overrightarrow{{F}_{i}}+\overrightarrow{{u}_{i}}\cdot \delta \overrightarrow{{F}_{i}}+\overrightarrow{{u}_{j}}\cdot \delta \overrightarrow{{F}_{i}}+\delta \overrightarrow{{u}_{i}}\cdot \delta \overrightarrow{{F}_{i}}+\delta \overrightarrow{{u}_{j}}\cdot \delta \overrightarrow{{F}_{i}}]{{\rm{d}}{S}}_{i}.\end{eqnarray} \tag{ 30 }$$

Here $\delta \overrightarrow{u}$ and $\delta \overrightarrow{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 $\underline{\underline{\tau }}$ that arises in the elastic environment of each force dipole in response to the displacement $\overrightarrow{u}$ on its surface. The energy evaluated in equation (30) equals the work done by spherical force dipole i on its surface, and thus we add the index i and the total interaction energy is ${\rm{\Delta }}{E}_{1}+{\rm{\Delta }}{E}_{2}$. The forces applied by each force dipole are equal and opposite to the forces applied on it by the environment: $\overrightarrow{F}=-\underline{\underline{\tau }}\cdot \hat{r}$, where $\hat{r}$ is the outward pointing unit vector normal to the surface of the sphere after its deformation and movement. Since we allow only volume change and translational motion of the entire force dipoles, the direction of the normal does not change due to the interaction. The terms of the stress tensor $\underline{\underline{\tau }}$ based on the general solution of the displacement field equations (5)–(6) [48] are given in appendix C. $\overrightarrow{{F}_{0}}=\tfrac{4{{Gu}}_{0}}{{R}_{0}}\hat{r}$ is the force per unit area on the surface of the sphere in the case of a single force dipole with known isotropic displacement $\overrightarrow{{u}_{0}}={u}_{0}\hat{r}$ on its surface and without interactions with other spheres, and thus the self energy is ${E}_{0}=8\pi {{Gu}}_{0}^{2}{R}_{0}$.

Evaluation of the integral in equation (28) may be simplified by noting that both the displacements $\overrightarrow{u}$ and the stresses $\underline{\underline{\tau }}$ are expressed in terms of Legendre polynomials and their derivatives (see equations (21)–(22)), and using their orthogonality we get (see appendix D):

$$\begin{eqnarray}&&E=2\cdot \displaystyle \frac{1}{2}{\int }_{S}\overrightarrow{u}\cdot \overrightarrow{F}{\rm{d}}{s}=2{E}_{0}\left(1+\displaystyle \frac{1}{\sqrt{4\pi }}{D}_{0}-\displaystyle \frac{1}{\sqrt{4\pi }}\left\{-{D}_{0}+\displaystyle \sum _{m=0}^{\infty }[{f}_{0m}^{{Cr}}{C}_{m}+{f}_{0m}^{{Dr}}({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})]\right\}\right).\end{eqnarray} \tag{ 31 }$$

Equation (31) is valid only asymptotically for $\widetilde{d}\gg 1$. The exact expressions for arbitrary distance and for all four regulation scenarios is given in appendix D. In the FS cases D₀ is evaluated from the cancellation condition and equation (31) becomes:

$$\begin{eqnarray}&&{E}_{{FS}}=2{E}_{0}\left(1+\displaystyle \frac{1}{\sqrt{4\pi }}{D}_{0}\right)=2{E}_{0}\left(1+\displaystyle \frac{1}{\sqrt{4\pi }}\left\{\displaystyle \sum _{m=0}^{\infty }[{f}_{0m}^{{Cr}}{C}_{m}+{f}_{0m}^{{Dr}}({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})]\right\}\right),\end{eqnarray} \tag{ 32 }$$

while in the VS cases D₀ is zero and equation (31) becomes:

$$\begin{eqnarray}&&{E}_{{VS}}=2{E}_{0}\left\{1-\displaystyle \frac{1}{\sqrt{4\pi }}\displaystyle \sum _{m=0}^{\infty }[{f}_{0m}^{{Cr}}{C}_{m}+{f}_{0m}^{{Dr}}({D}_{m}-\sqrt{4\pi }{\delta }_{m,0})]\right\}.\end{eqnarray} \tag{ 33 }$$

The meaning of the sum of the last two terms in equations (32)–(33) is the change of the volume of force dipole 1 caused by all the modes of displacement created by force dipole 2 (see figure 2), ${\rm{\Delta }}{V}_{21}=g$, see equation (26). Thus the interaction energies in the FS cases become:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{FS}}=2{E}_{0}{\rm{\Delta }}{V}_{21},\end{eqnarray} \tag{ 34 }$$

and in the VS cases:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{VS}}=-2{E}_{0}{\rm{\Delta }}{V}_{21}.\end{eqnarray} \tag{ 35 }$$

Note that in the VS cases ${\rm{\Delta }}{V}_{21}$ is the change in the symmetric (n = 0 mode) radial displacement due to the interaction, while in the FS cases ${\rm{\Delta }}{V}_{21}$ is the correction to that mode that maintains the total displacement unaffected by the interaction.