Many-body interactions between contracting living cells

Abstract The organization of live cells into tissues and their subsequent biological function involves inter-cell mechanical interactions, which are mediated by their elastic environment. To model this interaction, we consider cells as spherical active force dipoles surrounded by an unbounded elastic matrix. Even though we assume that this elastic medium responds linearly, each cell’s regulation of its mechanical activity leads to nonlinearities in the emergent interactions between cells. We study the many-body nature of these interactions by considering several geometries that include three or more cells. We show that for different regulatory behaviors of the cells’ activity, the total elastic energy stored in the medium differs from the superposition of all two-body interactions between pairs of cells within the system. Specifically, we find that the many-body interaction energy between cells that regulate their position is smaller than the sum of interactions between all pairs of cells in the system, while for cells that do not regulate their position, the many-body interaction is larger than the superposition prediction. Thus, such higher-order interactions should be considered when studying the mechanics of multiple cells in proximity. Graphic Abstract


Introduction
Live cells exert contractile forces on their environment.The shape, size, and resulting biological change in response to changes in the rigidity of the ECM [3][4][5].It is not fully clear how cells respond to changes in the mechanical environment caused by other cells, external forces, or changes in the rigidity of the medium.In many studies, the working hypothesis has been that cells tend to maintain specific quantities through mechanical homeostasis [6,7].For example, by regulating the forces they apply, cells will vary the displacements they generate as their environment changes.Alternatively, cells may change the forces needed to create those displacements by regulating their deformation.Furthermore, cells modulate their shape and spatial contractility patterns in response to environmental changes.
The mechanical activity of cells is often described by force dipoles, namely pairs of equal and opposite active forces that each cell applies on its mechanical environment [8][9][10].There are analogies between such force dipoles and electric dipoles that consist of two equal and opposite electric charges.Similarly, mechanical interactions between cells result from each cell generating a deformation field in the surrounding medium, which resembles the electric field formed around an electric dipole.Distant cells are, in turn, influenced by this field.Thus, matrix-mediated interactions between cells are similar but not identical to interactions between electric dipoles.
Fig. 1: Each cell is modeled as a spherical force dipole, comprised of radial active forces that are isotropically distributed on the surface of a sphere Nonetheless, one can study the elastic interaction between spherical cells surrounded by a linearly elastic material [7,13,23].The concepts introduced and the physical mechanisms identified in such studies are also relevant to morphologically complex cells in nonlinear materials.Specifically, despite the linear properties assumed for the ECM, the intra-cellular mechanisms for regulating each cell's mechanical activity give rise to nonlinearities that show up in inter-cellular behavior.
In this paper, we investigate how cellular regulation of mechanical activity breaks the superposition that one could naively expect to find due to the linear elastic response of the surrounding medium.We analyze situations containing multiple contractile cells and show that the total interacting energy in such cases differs from the result obtained by assuming that the interactions are pairwise additive.

Shape Regulation
We distinguish between two types of spherical force dipoles, based on the presence or absence of regulation of the forces that they apply; dead, but active force dipoles do not regulate the forces that they apply, and their activity does not depend, for instance on the distances to their neighbors.In our model, live cells are capable of measuring external forces and deformations on their surface and adjusting the active forces that they apply according to some internal algorithm, for example, to maintain a certain displacement or a certain force on their surface.
The difference between these two types of behavior is evident when a spherical cell generates a radial and isotropic self-displacement field, i.e., the displacements induced by this cell in the absence of neighboring cells.Such a field would cause the cell to only change its volume, without any distortion of its shape.In that case, a pair of such dead active force dipoles preserves their selfdisplacement fields, leading to vanishing interaction energy [23].Note that if the self-displacement field generated by each sphere is radially symmetric around the center of that sphere, then the total displacement on the surface of each cell, which is the sum of the self-displacement fields generated by these two contracting spheres, would be anisotropic, as shown in Fig 2a .This result is general for contracting objects or arbitrary shapes, that generate self-displacements only in their principal directions [24,25].(b-c) Live force dipoles regulating the force they apply to remain spherical even in the presence of the other force dipole.Initial shape (purple), forces applied by each force dipole (red arrows), corresponding self displacements (dashed black), displacements caused by the other force dipole (dashed green), total displacement (blue), the center of the interaction-free displacement (black dot), the center of the displacements with the interactions (blue dot).For illustration purposes, the initial distance between the spheres was set to d = 3R 0 , and the self-displacement to u 0 = 0.4R 0 , where R 0 is the radius of each sphere.The Poisson ratio is ν = 0.3.See text for description of the different regulation scenarios in (b-c) In comparison, two live cells that adjust the self-displacement fields that they generate have non-vanishing interaction energy.We focus on live cells with shape regulation [23], as demonstrated in Fig. 2b,c.Namely, spherical cells that adjust the anisotropic azimuthal distribution of the active forces that they apply, such that the total displacement on their surface will be radially symmetric.This total displacement is the sum of the selfdisplacement that each cell generates on its surface plus the displacement fields on its surface due to the activity of the other cells in the system.

Many-body Interactions
We consider a series of identical live spherical cells of radius R 0 , arranged along a straight line, separated by equal distances d between their centers, and surrounded by a three-dimensional linear elastic material of bulk modulus K and shear modulus G.In the absence of other cells, each cell contracts isotropically with a displacement u 0 on its surface.We assume that each cell senses the displacements created on its surface by all other cells and adjusts its active force to compensate for them and prevent distortions of its spherical shape.For this complex situation of multiple activity-regulating cells, we will calculate the total elastic energy stored in the surrounding medium.We will then obtain the many-body interaction energy by subtracting from the energy of this mutual situation, the sum of the self-energies of the cells, namely the energy stored in the medium assuming that each cell contracts independently in an infinite medium, without any other cells.As we will show below, due to the nonlinearity in the regulation of the active forces that the cells apply, this manybody interaction energy differs from the sum of all two-cell interactions.We develop a formalism for an arbitrary number of cells, and will explicitly solve the geometries including three or four cells, and compare them to the pair-wise additive result obtained from the analysis of geometries of two cells.We will also consider an infinite array of equally spaced cells along a straight line, for which we will calculate the interaction energy per cell.

Displacements Created by Spherical Cells
We describe the displacements generated by cells as the sum of an isotropic constant displacement u 0 and an anisotropic, interaction-dependent displacement ∆u that is intended to cancel anisotropic displacements caused by other cells.
The symmetry of the arrangement dictates that the displacement fields produced by cells placed at equal distances on either side of the array's center are mirror images.For simplicity, we place the origin of the coordinate system at this center and number the cells according to their distance from it, see Fig. 3.We choose the coordinate systems of the cells based on their index: left-handed for positive indices and right-handed for negative or zero indices, see Figs. 3, 4.This choice of coordinate systems is based on the system's symmetry and will simplify the calculations.
Fig. 4: Three spherical cells each with radius R 0 , all applying a radial isotropic displacement u 0 (red arrows) on their surfaces.The coordinate systems of spheres 0 and −1 are right-handed (blue and green accordingly) and the coordinate system of sphere 1 is left-handed (magenta) and may be written as θ 1 = π − θ ′ 1 where θ ′ 1 is the commonly-used right-handed azimuthal coordinate for sphere 1 The displacement field u around each cell must satisfy mechanical equilibrium [26]: Due to the rotational symmetry around the line passing through the centers of the cells, 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 cylindricallysymmetric displacements on its surface [26], we write the anisotropic displacements field satisfying Eqs.(2-3) outside the cell (r > R 0 ) as a multipole expansion in terms of spherical harmonics Y n (θ) = 2n+1 4π P n (cos θ): with the Legendre polynomial of order n [27].Equations (5-6) represent the anisotropic displacements created by each cell in its coordinate system with its origin in its center.Here, u ri and u θi are the radial and angular components of the displacement field caused by cell i, and the infinite sums represent the anisotropic corrections that each cell produces to cancel the shape distortion caused by its neighbors.We have inserted u 0 and R 0 to make the coefficients C i n and D i n dimensionless.Using the dimensionless displacements u ri = uri u0 , u θi = u θi u0 , and the dimensionless position r = r R0 , we rewrite Eqs. ( 5) and ( 6) as follows: Note that Eqs.(8-9) solve Eq. ( 1) only when each cell is surrounded by an infinite, homogeneous linearly-elastic medium, including in the interior of the neighboring cells.Biological cells have a rigidity that differs from the rigidity of the ECM that surrounds them; thus, this assumption seems 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 of how the cell generates this stress on its surface.In particular, the stress that actual cells apply on their surrounding includes passive stress coming from the rigidity of the cell plus 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 it performs, which determines the interaction energy, and our results are valid irrespective of the mechanical rigidity of the cells themselves.See also Ref. [13].

Cancellation Condition
To preserve isotropic displacements on their surface, live cells in our model create correcting displacements that cancel the anisotropic displacements created by their neighbors.Thus, the sum of all anisotropic displacements caused at the surface of a cell by all other cells and all the corrections applied by the discussed cell must vanish.The coefficients C n and D n in Eqs.(5)(6) are derived in this way so that each cell can retain its spherical shape despite interacting with its neighbors.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 fields of each cell j to the coordinate system of the discussed cell i by substitution of the expressions for r j and θ j in terms of r i and θ i and then multiplying the displacement vector − → u j = (u rj , u θj ) by a transformation matrix.The transformation matrix depends on the coordinate systems of the cells i and j; we use the rotation matrix for i and j with the same signs, and the reflection matrix for indices with opposite signs.The central cell i = 0 may be treated as having a positive or a negative sign and right or left-handed coordinate system, accordingly.In our analysis, we chose to treat the central cell as having a left-handed coordinate system, and thus, we treat its index as positive (see Fig. 3).
We write the resultant expressions for the radial and angular displacements caused by each cell j on the surface of cell i in terms of the spherical harmonics of cell i by writing the projections: As may be seen from Eqs. (12,13), every sphericalharmonic mode of cell j contributes to all the modes on the surface of cell i.Finally, we sum the contribution from all cells and find the total displacement in each mode.The anisotropic displacements caused by all cells j = i must be canceled on the surface of cell i by the corrections it applies.We write the dimensionless displacement u ii created by cell i on its surface (namely at r i = 1): The term δ n,0 in Eq. ( 14) is a Kronecker delta, which represents the isotropic radial displacement created by cell i on its surface without the anisotropic cancellation corrections.This constant term does not depend on changes in the cell's environment.The remaining terms are different modes of additional displacement that this cell creates in response to the displacement field induced on its surface by the neighboring cells.The dimensionless displacement u ji created by each cell j on the surface of cell i is: where the sum over m originates from the fact that the displacement u j created by cell j is given by a multipole expansion (8-9) with the corrective magnitudes C j m and D j 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 cell i, each mode from cell j contributes to all the modes of cell i. depend only on the dimensionless distance d ji = dji R0 between the cells.However, similarly to the transformation matrices, these functions depend on whether the indices i and j have the same or opposite signs.This follows from the choice of the coordinate systems of the cells that were described earlier, see Fig. 3.If the signs are the same, the functions further depend on the sign of the difference |i|−|j| that indicates the side at which cell j is located relative to cell i.Thus we make a distinction between f Cr nm l , f Dr nm l , f Cθ nm l , f Dθ nm l and f Cr nm r , f Dr nm r , f Cθ nm r , f Dθ nm r for |j|> |i| and |j|< |i|, accordingly, with i and j of the same signs, and f Cr nm o , f Dr nm o , f Cθ nm o , f Dθ nm o for i and j with opposite signs.The expressions for all cases are given in Appendix A.
We now require that for live cells, the total displacement u ii (θ i ) + j u ji (θ i ) on the surface of cell i 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 cells.Moreover, we require that the center of symmetry of each cell does not move.This will be denoted fixed size fixed position (FSFP) regulation.We will also consider three additional activity regulation scenarios in which the interaction causes the cells to change their volume and/or to move, yet they remain spherically symmetric.We denote these regulation scenarios as: variable size fixed position (VSFP), fixed size variable position (FSVP), and variable size variable position (VSVP), see Fig. 5 and Fig. 2 above.For VSFP, cells regulate their shape and rigid body motion, but not their size.In this case, we nullify D 0 , the first term of the regulating series of each cell responsible for cell size regulation.Similarly, for FSVP, in which cells regulate their shape and size but not their position, we nullify the coefficients C 1 and D 1 , and for VSVP, in which cells regulate their shapes but not their size or position, we nullify D 0 , C 1 , and D 1 .This method is discussed in further detail in Ref. [23].
To preserve isotropic displacement on the surface of cell i we require that: Due to the symmetry of the system and our choice of coordinate systems for the cells, the coefficients of pairs of cells with opposite indices are equal, namely C j m = C −j m and D j m = D −j m .Thus, for a system of k cells, we need to write the conditions (18,19) for k/2 cells with a nonrepeating index j if the total number of the cells is even, and for k/2 + 1 if it is odd.
Substituting (14,15,16,17) in (18,19) yields: Fixed Size: Variable Size: Fixed Position: Variable Position: h g g h Fig. 5: Schematic drawings showing how the four possible scenarios of shape regulation depend on whether the size (mode n = 0) and the position (mode n = 1) are regulated or not.Initial shape (solid purple), displacement without (dashed black), and with (solid blue) interaction.In the variable size cases, the cell is maintained spherical, yet its radius changes by g.In the variable position cases, the center of the cell (dot) translates by h 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 ≥ 1 we require: Note that for n = 0, from Eqs. (8-9) C 0 is irrelevant; thus, we set it to zero.Moreover, since = 0 and Eq. ( 21) holds trivially, thus for n = 0 we obtain only one equation, from Eq. (20): We obtain closure of the infinite coupled linear Eqs.(22)(23)(24) by assuming that C n = 0 and D n = 0 for n > 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 cells, and since the solutions decay as 1/r n , at large r, large n terms become negligible.We previously verified this numerically by increasing n max until convergence [23].According to our findings, for two cells, n max = 1 for FP regulation and n max = 2 for VP regulation scenarios are sufficient to include the leading terms and to obtain good approximations for the interaction energy.Therefore, we use these values of n max also in the present analysis of interactions between multiple cells.The resultant expressions for the coefficients C n and D n for three and four cells along a straight line are given in Appendix B. We evaluate the forces created by the cell using Eq.(C22-C27) in Appendix C for the stress tensor in the elastic environment.Due to force balance, the cell's active force per unit area is exactly minus this elastic stress.
The case of many cells along a straight line can be solved by approximating it by an infinite, one-dimensional array of cells; an infinite number of neighbors surrounds each cell.Consequently, all cells respond similarly to their environment and create identical displacement fields.Therefore, C i n = C j n and D i n = D j n for any i and j.This reduces the number of unknown coefficients from n max × k in a finite array of k cells to n max in an infinite array, enabling us to define cancellation conditions for a single general cell rather than for k/2 cells.
To be finite and solvable, we include in Eqs.(18,19) only terms coming from a limited number k max of neighboring cells, despite the assumption that there is an infinite number of cells.Similarly to n max , this is justified since the displacement fields created by the cells decay as 1/r n , so displacements produced by distant cells become negligible.As shown below, we verify this numerically by increasing k max until convergence.

Interaction Energy
We solved the equations for configurations of two, three, four, and an infinite number of cells on a straight line.For each case, we evaluated the extra work performed by each cell i by terminating the infinite sums at n max = 1 for FP regulation and at n max = 2 for VP regulation.For configurations that involved three or more cells, we compared the interaction energy obtained from the direct solution of the multiple-cell geometry with the pair-wise additive prediction assuming superposition of interactions between all pairs of cells within the system.The direct calculation consists of constructing and solving a set of boundary conditions of the form of Eqs.(20)(21).The superposition calculation approximates three-, four-, and manybody interactions by summing all the two-cell interactions in the system.
The total elastic energy stored in the medium surrounding the cells is equal to the work performed by all cells to generate their deformations, starting from their undeformed states.Cells apply active forces only on their surfaces, thus the amount of work performed by each cell at any point on its surface may be computed by multiplying the force that the cell applies at that point by the total displacement there, divided by two.The division by two results from the integration starting from the undeformed state and reaching the deformed state as the stress in the system gradually builds up linearly with the growing displacement in our linearly elastic medium [13].The self-energy of each cell is the elastic energy it generates when it is surrounded by the infinite ECM and is isolated from other cells.We define the interaction energy as the difference between the elastic energy of the system of cells and the sum of all the cells' self-energies and is thus equal to the extra work performed by the cells due to the presence of other cells around them.
We write the extra work performed by cell i in a case that includes k cells as: where E 0 = 8πGu 2 0 R 0 is the cell's self-energy, or the work done by a single, isolated cell that creates on its surface an isotropic displacement u 0 [23].We find that for all the cases that we considered, the dimensionless extra work may be written as Here, A k i is a numerical prefactor that depends on the number k of cells in the system and on the index i of the cell within the system, but which does not depend on the medium's Poisson ratio ν.We find that this dependence may be included in B(ν), which is the same for all cells within the system and for any number of cells in the system.Finally, α is the exponent of the power law decay of the interaction energy with the distance between the cells.We find that the sign of the extra work is q = +1 for FS and q = −1 for VS.These signs are consistent with the theoretical understanding that FS refers to displacement homeostasis, which leads to repulsion between cells and VS to stress homeostasis, which leads to attraction [7,13,23].Table 1 shows the values of A k i , α, and B for different cells in configurations that include different numbers of cells on a straight line, and for the different position regulation scenarios.Note that by symmetry, A k j = A k −j .For an infinite array of cells, all cells are equivalent.This symmetry cancels the displacements produced by the cell's neighbors, thus its position remains fixed even without position regulation, and position-regulating terms vanish in all regulation scenarios.Thus, the distinction between FP and VP becomes irrelevant.Nonetheless, we refer to the results here as VP regulation, since the positions of the cells are not actively regulated by the cells similar to the VP scenarios in the two, three, and four cell configurations.
We summed these additional works per cell to evaluate the direct interaction energy for a configuration of k cells, where A k tot = i A k i is also given in Table 1.
We find that the sign of the interaction energy, given by q, as well as the scaling with distance, given by the exponent α, are independent of the number of cells.Moreover, for the three-and fourcell cases, the additional work performed by the central cells is greater than that performed by the side cells, namely A 3 0 > A 3 1 and A 4 1 > A 4 2 .This is due to the fact that the central cells are closer to the rest of the cells compared to the side cells.Using the direct method, we did not find a closed-form solution for an infinite array of cells on a straight line.Thus, the number of interacting neighbors included in computations, in this case, is limited.Nevertheless, the addition of interacting neighbors does not influence the sign q of the interaction or the scaling exponent α with cellcell distances.The energy remains proportional to d −6 , and only the coefficient A ∞ i is affected by the number of neighbors included in the calculation.Figure 6 shows how the value of A ∞ i converges as the number of interacting cells grows, and the value given in Table 1 is for the largest number of cells that we considered, k max = 39.If the interaction energy was pair-wise additive, one could treat the interactions of three, four, and many cells as combinations of two-cell interactions between all the cells in each configuration.The total interaction energy would then be equal to the sum of energies of interactions between all pairs of cells and may be evaluated using the two-cell results given in Table 1.For example, we consider the interaction energy between three cells in the FSFP case by decomposing it into two similar interactions between the side cells and the central cell and the interaction between the cells on opposite sides.Since the distance between these cells equals twice the distance d between the side cell and the central cell, the interaction energy becomes: In the same manner, we evaluate the added work performed by a cell in an infinite array of cells as a sum of pair interactions with each cell on both sides, where according to Table 1, in the VP cases, the interaction energy of each pair equals We denote the results obtained from this superposition calculation by the subscript s, and list these results as well in Table 1.We denote the difference between the direct many-body calculation and the superposition expression as ∆A ≡ A tot − A s .
From Table 1, we see that the coefficients A i found using the two methods for the same configurations are different.This difference is not obvious since, in linear elasticity, typically results may be superimposed when analyzing more complicated arrangements.We conclude that the active response of the cells to their neighbors produces non-linear intercellular interactions, even in the case of linear elasticity.Considering that each active cell in the presence of other active cells is performing extra work, one might expect the interaction energy in all cases to be higher in direct method solutions than in superposition method solutions.However, in three-and four-cell configurations, this assumption is correct in the VP but not in the FP regulation scenarios.Namely, for VP, A tot > A s , while for FP A tot < A s .This unexpected result follows from canceling the central cell's rigid body motion due to the configuration's symmetry.For FP, a large part of the added work comes from the interaction between the size regulation of the cell (n = 0 mode) and the forces that regulate the motion of the neighbor (n = 1 mode) as a rigid body.The added work done by the central cell is significant due to its interaction with two neighbors on both sides.The work done by the side cells is small due to the absence of firstmode forces created by the central cell, see Fig. 7a.Most of the added work done by the side cells follows from their interaction with the relatively distant cells on the opposite sides and is small due to the fourth power of normalized distance d.In  For FP, the cells on the sides apply forces to regulate their motion.Due to the symmetry, the central cell does not need to apply forces to stay in place.For VP, no cell regulates its motion.Initial shape (purple), forces applied by each cell (red arrows), corresponding self displacements (dashed black), displacements caused by the other cells (dashed green), total displacement (blue), the center of the interaction-free displacement (black dot), the center of the displacements with the interactions (blue dot).For illustration purposes, the initial distance between the spheres was set to d = 3R 0 , and the self-displacement to u 0 = 0.4R 0 .The Poisson ratio is ν = 0.3 contrast to FP cases, position-regulating terms are assumed to vanish and are not affected by symmetry in VP cases, see Fig. 7b.The first mode is the only one that includes an antisymmetric function; thus, only this mode will be affected by symmetry.

Conclusions
We model live cells in the ECM as spherical active force dipoles, which are surrounded by a linear elastic environment.For isotropic active forces and thus isotropic self-displacements, the interaction energy between cells vanishes.Hence, we distinguish between this dead behavior, in which the cells apply constant forces and self-displacements on their surface, and live, regulatory behavior, in which cells adjust their active forces and selfdisplacements in response to changes that they sense in their environment.This live behavior of cells is similar to interaction between induced electric dipoles on particles with charge regulation.
We examine systems with three, four, and infinite numbers of cells on a line.We solved the interaction energy for these configurations for four different types of self-regulation: on top of preserving their spherical shape, cells can also preserve their volume or their position, or both.Similarly to the interaction between two such shape-regulating cells [23], for fixed position, we found the interaction energy to be inversely proportional to the distance between the cells to the fourth power, and for variable position, to its sixth power.As in the case of two cells, also here, we found that for fixed volume, multiple cells are repelled from each other, and for variable volume they are attracted to each other.
We compared the results of direct computation of the many-body configurations to the sum of all two-cell interactions for the same configurations.A comparison of the results shows that the superposition method does not predict the energy of multiple-cell configurations.We also found that if cells regulate their position, the many-body interaction energy is smaller than the sum of interactions between all pairs of cells in the system, while for cells that do not regulate their position, the many-body interaction is larger than the superposition prediction.We conclude that the active response of cells to their neighbors produces non-linear intercellular connections even in the case linear elasticity.
We have solved the deformation fields for the case in which the rigidity of the cells is the same as that of their environment.Biological cells, however, are complex entities whose rigidity varies from place to place and from the rigidity of the ECM.To relate our results to live cells, we describe each of them as a mechanism that applies forces on the surface and responds by their variation to the application of external force or displacement.The displacements and forces applied by a cell may be divided into "dead" and "live" parts.While the dead part of the forces or displacements would remain the same if the cells were dead and retained their elastic properties, the live part depends on their programmed behavior and is generated by the contraction of their actomyosin networks.Since the resultant force and displacement are the sum of those two parts, cells may create such a live response so that the resulting forces and displacements will coincide with the case considered here, for which their rigidity is identical to the rigidity of their environment.Even if cells do not behave in this manner, our results highlight the many-body nature of matrix-mediated elastic interactions between cells, and specifically the different behavior for different regulation scenarios.
Following our work on multiple cells along a straight line, it would be interesting to extend our work to two-dimensional arrangements, the simplest of which would be three cells at the corners of a triangle.Since there will be no cylindrical symmetry like in our present study, a more complicated analysis of displacements on the surface of the cells will be needed to model their behavior.It would also be interesting to expand our present work to the case of aspherical cells, for example, oblate spheroids.In this case, the interaction energy between two such cells would depend on the distance between their centers and on the relative angle between their axes.We limited ourselves to cells surrounded by a linearly elastic medium, so that we could exactly solve their interactions analytically.It would be interesting to test our qualitative predictions by solving with numerical simulations situations with nonlinear response of the medium.
Appendix A The functions f Cr nm t , f Dr nm t , f Cθ nm t and f Dθ nm t We distinct between three different cases for the functions f Cr nm t , f Dr nm t , f Cθ nm t and f Dθ nm t appearing in Eqs.(16,17).Thus the index t in the following expressions may be equal to l, r or o: We used the identity [28]: to rewrite the derivatives dYn(θ) dθ in terms of θ, and for the sake of brevity we defined:

Appendix C The stress tensor
Here, we develop the expressions for the stress tensor in the cases discussed in the paper.The applied displacements are symmetric about an axis passing through the centers of the cells.The expressions are taken from [26] for the case of the displacement field given by Eqs.(5-6), excluding the first term in Eq. ( 5), which corresponds to volume change.The stress tensor may be written in the following form in this case: The stress tensor is symmetric τ ij = τ ji and thus only six components are to be evaluated.We define the dimensionless stress tensor τ = τ G R0 u0 , the elements of which are given by: τ   After obtaining the coefficients C n and D n , we compute the extra work W k i performed by cell i in a configuration that consists of k cells, to generate total displacement − → u in accordance with Eqs.(14)(15)(16)(17) and the type of the regulation: Here, the integration is over the spherical surface of the cell i, − → F is the force per unit area applied by it on its environment.Due to force balance, the active forces applied by each cell are equal and opposite to the forces applied on it by the environment: − → F = −τ • r, where r is the outward pointing unit vector normal to the surface of the cell after its deformation and movement, and τ is the stress tensor given above that arises in the elastic environment of each live cell in response to the total displacement − → u on its surface.Due to the spherical shape regulation of the cells, there are no displacements in the azimuthal direction in this study.
− → F 0 = 4Gu0 R0 r is the force per unit area on the surface of a single cell with known isotropic displacement − → u 0 = u 0 r on its surface and without interactions with other spheres, and thus − → u 0 • − → F 0 = 8πGu 2 0 R 0 .

Fig. 2 :
Fig.2: Two spherical force dipoles: (a) Dead force dipoles applying an isotropic elastic force, (b-c) Live force dipoles regulating the force they apply to remain spherical even in the presence of the other force dipole.Initial shape (purple), forces applied by each force dipole (red arrows), corresponding self displacements (dashed black), displacements caused by the other force dipole (dashed green), total displacement (blue), the center of the interaction-free displacement (black dot), the center of the displacements with the interactions (blue dot).For illustration purposes, the initial distance between the spheres was set to d = 3R 0 , and the self-displacement to u 0 = 0.4R 0 , where R 0 is the radius of each sphere.The Poisson ratio is ν = 0.3.See text for description of the different regulation scenarios in (b-c)

Fig. 3 :
Fig. 3: Left-handed coordinate systems chosen for cells with positive index n > 0 and right-handed for zero or negative index n ≤ 0 for (a) odd number of spheres, (b) even number of spheres

Fig. 6 :
Fig. 6: Convergence of the coefficient A ∞i with the number of cells k max included in the calculation for an infinite array of cells.Exact evaluation (red squares) taking into account interactions of 3 to 39 cells.Approximate fit (dashed line) given by A ∞ i = 84.42x 4 +28.88x 3 −99.56x 2 +0.1618x+28.90where x = 1/k max .
(a) Interaction between cells in FP regulation case.(b) Interaction between cells in VP regulation case.

Fig. 7 :
Fig.7: Interactions between three cells in FP (a) and VP (b) regulation cases.For FP, the cells on the sides apply forces to regulate their motion.Due to the symmetry, the central cell does not need to apply forces to stay in place.For VP, no cell regulates its motion.Initial shape (purple), forces applied by each cell (red arrows), corresponding self displacements (dashed black), displacements caused by the other cells (dashed green), total displacement (blue), the center of the interaction-free displacement (black dot), the center of the displacements with the interactions (blue dot).For illustration purposes, the initial distance between the spheres was set to d = 3R 0 , and the self-displacement to u 0 = 0.4R 0 .The Poisson ratio is ν = 0.3 n (side cells) and C 0 n , D 0 n (central cell) of the multipole expansion at asymptotically long distances, d ≫ 1 in the case of interaction of three cells.Here Λ 1 = 1−2ν 5−6ν , Λ 2 = 1 5−6ν , Λ 3 = 1 4−5ν , n max is the highest-order term taken into account, and d = d R0 is the dimensionless distance between neighboring cells.

Table 1 :
Coefficients for additional work Eq.(26) done by the cells as a result of mechanical interaction with their neighbors, for different scenarios of position regulation.