Mechanobiological induction of long-range contractility by diffusing biomolecules and size scaling in cell assemblies

Mechanobiological studies of cell assemblies have generally focused on cells that are, in principle, identical. Here we predict theoretically the effect on cells in culture of locally introduced biochemical signals that diffuse and locally induce cytoskeletal contractility which is initially small. In steady-state, both the concentration profile of the signaling molecule as well as the contractility profile of the cell assembly are inhomogeneous, with a characteristic length that can be of the order of the system size. The long-range nature of this state originates in the elastic interactions of contractile cells (similar to long-range “macroscopic modes” in non-living elastic inclusions) and the non-linear diffusion of the signaling molecules, here termed mechanogens. We suggest model experiments on cell assemblies on substrates that can test the theory as a prelude to its applicability in embryo development where spatial gradients of morphogens initiate cellular development.

where V (x, x ) (which has units of inverse volume) accounts for the generally, long-ranged mechanical interactions of the contractile cells so that a cell at a position x affects a cell at position x arbitrarily far away through the strain it induces in the underlying elastic medium. This accounts for the possibility that the local contractility is also actively changed by the cell in response to the local strain induced by the other contractile cells, such as through mechanically transduced myosin regulation 2 . This linear response of the local cellular contractility to both the mechanogen concentration and the contractility of the rest of the system can be motivated by considering the mechanical energy of the combined system of contractile cells and the mechanogens. We include in the energy terms the interaction of the cellular force dipoles 1 with each other as well as the binding of mechanogens to them as, where ψ(x) represents the local contractility (trace of the force dipole density tensor) and c(x, t) the concentration of mechanogens; the energy is scaled by the energy cost of contractility in the absence of either other cells or mechanogens (the coefficient of the first term in Eq. (2)), and F c [ψ(x), c(x, t)] = −χ dx c(x, t)ψ(x) represents the tendency of mechanogens to induce contractility with χ > 0. Since mechanical forces are transmitted very quickly 3 , the cell assembly can be considered to be in mechanical equilibrium and local binding equilibrium even if the mechanogens are not yet in their global, steady-state concentration profile. The contractility depends in an adiabatic manner on the local mechanogen concentration at a given time, c(x, t), which evolves much more slowly than the dynamics of the mechanics. Therefore, minimizing the mechanical energy in Eq. (2) with respect to the contractility ψ, yields the condition for mechanical equilibrium, Eq. (1). However, this is only an example that motivates the linear response relation, Eq. (1), which is more general. We close the mechanogen-contractility feedback loop by considering a scenario where the degradation of mechanogens is due to the cell contractility, (i.e., the trace of the dipole tensor) which can happen through changes in their rate of uptake (via changes in membrane tension) or through cooperative binding of the mechanogens that show a preference to bind to contractile cells: This corresponds to case (ii)of strain-determined degradation considered in the main text, where the degradation couples to the local contractility instead of strain. Eqs. (1) and (3) together then provide a different, more formal, formulation of the mechanogen-contractility coupling that does not make any explicit reference to the strain induced by the dipoles, and are an alternative to the model described in the main text in terms of Eqs. (1) and (2). The effective equation we get for the mechanogen concentration after integrating out the contractility, ψ(x) is equivalent to Eq. (4) in the main text, except for a reversal in the dependence on the mechanical boundary conditions, or in other words a change in the sign of the linear degradation term in Eq. (4). As we show in the main text, the scaling relations obtained for the decay length of the mechanogen profile with system size are not affected by this. Thus the size scaling is quite robust to the exact model or mechanism used for coupling mechanogen concentration to cell mechanics.

II. INTERACTIONS OF FORCE DIPOLES IN A FINITE MEDIUM
The theory of continuous elastic media determines that a distribution of forces induces stresses in the medium that satisfy a force balance 4 as, where σ ij is the stress tensor (force per unit area), and f i (x) the local density of forces (per unit volume). In a linear elastic medium, the stress depends linearly on the strain through the elastic tensor, where u kl (x) = (1/2)(∂ k u l (x) + ∂ l u k (x)) is the strain tensor in terms of the local displacement u(x) of the medium. We now relate the local force dipole density p ij (x) (dipole moment per unit volume) to the force density f (x). Force dipoles comprise two equal and opposite forces F that are separated by a distance a. The product of the force and the distance defines a local force dipole 3 In a continuum representation (valid for scales much larger than |a|), one can consider the force density f (x) which is zero when integrated over some small volume ∆ due to the equal and opposite forces of the dipoles. The force dipole is the first moment of the force density integrated over that volume.
This relation can be inverted to give the local dipole density p ij (x). We consider a force distribution that comprises N pairs of equal and opposite forces separated by corresponding small distances (force dipoles). Writing this explicitly and expanding the expression for small distances, a between the forces we get where we chose the signs of the δ functions so anti-parallel alignment of the vectors a and F result in a contractile force dipole. The force density 8 can then be written as a force dipole density, where ← → p is the force dipole density tensor and the value of the local dipole tensor of the m th dipole is The work done by force dipoles to deform the medium is given by the product of the force applied by the dipole and the local displacement of the medium that is caused by that force, summed over the whole volume, as which using Eq. (11) can be written in terms of the force dipole density as, This after integrating by parts can be expressed as, which has contributions from the dipole density in the bulk as well as from the surface of the finite domain, wheren is the unit surface normal. This is expected as a distribution of force dipoles induces both volume and surface forces.
(Here V implies an integral over the volume and S an area integral over the surface of the finite medium.) The boundaries could be clamped or held fixed, such that the displacement is constrained to be zero at the surface, thus implying that the surface integral in Eq. (14) above vanishes. The other physically relevant condition is when the boundaries are free to move, or in other words are stress-free. This can be realized only if the surface forces from the dipoles in the medium are balanced by an externally imposed surface force distribution from the region outside the finite domain under consideration. Because the surface terms are zero in both of these two cases (Dirichlet and Neumann), the work done by a distribution of dipoles on the medium is the product of the dipole density with the local strain field 1,6 , The central point of this Appendix is that this energy can be interpreted as the interaction of a dipole at a position x with the local strain field which is induced by all the other dipoles in the system 1,6 .

III. ELASTIC DIPOLE INTERACTIONS IN A FINITE 1D DOMAIN
We now derive the mechanical interaction energy of a pair of force dipoles in 1D by considering an intuitive picture of a line of springs connected in series. A force dipole is represented as a single spring with a pair of equal and opposite forces (in this case compressive) applied on either end of it. We describe the discrete set of springs with a coarsegrained field u(x) which characterizes the displacement of the springs at a position x. The i th spring is considered to lie between x = (i − 1) · l and x = i · l where l is the length of each spring which is also the coarse-graining displacement in this description. We now apply this to a model of a finite medium, represented by a total of N springs thus forming a line of length L = N l. The mechanical equilibrium condition relates the change of the deformation with a general discrete distribution of forces as 4 , where for a spring system, the 1D elastic modulus E ≡ kl, is more conveniently expressed in terms of the spring constant k (units of force per length). A single force applied at a point x m therefore causes a jump discontinuity (seen by integrating Eq. (16)) in the slope of the displacement field, i.e. in the strain field, This just implies that the springs to one side (right or left) of the force are compressed while those on the other side are stretched. A force dipole that compresses the m th spring therefore implies two equal and opposite jumps in the slope of the displacement at positions (m − 1) · l and m · l which causes all the springs to the left or right of the "dipole spring" to be deformed equally. There is no deformation of the other "external springs" if the ends of the line of springs are free to move since in that case the springs just translate. However, a nontrivial deformation can arise if the line of springs is held fixed or clamped at either end such that we require u(0) = u(L) = 0. For a compressive force dipole that consists of a pair of forces, each of magnitude F , applied to either end of the m th spring, we find that each of the other springs in this line is stretched by a distance F/(N k) while the m th or "dipole" spring itself is compressed (change in length) by an amount F/k, thereby maintaining the constant overall length of the line. This solution satisfies the continuity of the displacement and the (two equal and opposite) jumps in its slope at (m − 1) · l and m · l respectively as suggested by Eq. (17). The deformation of the dipole spring is much higher than those of the other springs representing the elastic medium (by a factor of N ). In other words, if the dipole applies a pair of forces, each of magnitude F , the force in each of the other springs goes as F/N . The total energy stored in these springs then scales as N × (F/N ) 2 ∼ 1/N , i.e. inversely as the length of the line of springs. Let us consider next two force dipoles applied at the m th and n th springs where m and n are arbitrarily far apart. Let the forces of each dipole have magnitude F 1 and F 2 respectively. It can be shown by solving Eq. (16) with jumps in slope at positions m, m + 1, n, and n + 1 given by Eq. (17), that the deformation (change in length) of each of the "external" springs, i.e. all springs excluding the m th and n th ones, is (F 1 + F 2 )/(N k). Summing over the deformation energies of all these springs yields the total energy of the line of springs with two dipoles. Subtracting from this the sum of the total energy in the line of springs for a single dipole at m and at n respectively, we obtain the energy of interaction of the two dipoles as, which is independent of the relative location of m and n, is negative, and depends inversely on the length of the line of springs. This intuitive description in terms of discrete springs thus confirms the form of interaction of force dipoles suggested by Eq.(4) in the main text. Note that we have used the work done to deform the springs rather than the energy stored in them as our definition of energy since this is the more relevant quantity for an active dipole 1 .

IV. ELASTIC DIPOLE INTERACTIONS IN FINITE 3D SYSTEM
For an isotropic elastic medium, the condition for mechanical equilibrium for a distribution of forces can be written in terms of two elastic parameters as 4 (this is the higher dimensional analog of the 1D mechanical equilibrium expressed in Eq. (16)), where u(x) is the displacement of the medium at a position x induced by the local force density f (x), and K and µ are the usual bulk and shear modulus of the material. If the forces arise from an isotropic distribution of spherical force dipoles p ij (r) = p(r)δ ij , the forces produced by them and the consequent displacements in the medium are in the radial direction. In spherical polar coordinates, the dipole density tensor components are: p rr = p θθ = p φφ = p, and all the other (off-diagonal)components are zero. In other words, the force dipoles have a trace(hydrostatic) component 3p(r) but no shear (deviatoric) components. Writing Eq. (19) in spherical polar coordinates 4,5 , we obtain, where we have used the relation between force and dipole density, f (r) = −dp rr /dr = −dp/dr, obtained in Eq. (11). The left hand side of Eq. (20) can be expressed as, where the prime refers to differentiation with respect to r. Integrating this expression twice in succession, we obtain the displacement as, where the integral is indefinite, and c 1 and c 2 are integration constants determined by suitable boundary conditions. Since symmetry requires that the system is fixed at the center r = 0, c 1 is determined such that the radial displacement can be written as, while the corresponding expression for the radial component of the local strain is, We consider this solution in a sphere of finite size R (corresponding to a finite spherical assembly of cells). The constant c 2 is determined from the boundary conditions at the surface of the sphere. If the sphere is clamped at the surface, i.e., u(R) = 0, the corresponding displacement can be written as, while absence of radial stresses at the surface (free boundary conditions) leads to a displacement, The trace of the strain(in radial coordinates) 4 defined as T r(u) ≡ du/dr + 2u/r, for clamped boundaries is then, and for free boundaries, where p is the volume average of the isotropic dipole density, and for the free boundaries, we have dropped the constant p(R) which would anyway not contribute to the elastic interaction energy. Exactly analogous to the 1D case, the local strain in the medium thus has two contributions: one from the local dipole density and the other being a sum of the contributions of all the other dipoles in the finite medium. Note that we use ψ to represent the nondimensional isotropic cell contractility in the main text; which is related to the force dipole density, here denoted by p as ψ(x) = p(x)/p max p max is the maximum contractility of a cell and p is negative by convention for contractile dipoles. The local dipole is always contractile and thus contributes a compressive strain irrespective of the boundary conditions. The other dipoles stretch the medium when the boundaries are fixed while they compress the medium if the boundaries are free. This is reflected in the change in sign of the second term between Eqs. (27) and (28). The mechanical energy invested by force dipoles in deforming the medium was derived in Appendix II to be the product of the local dipole density and the local strain field produced by all the other dipoles, which is represented by the term −p/(K + (4/3)µ) in Eqs. (27) and Eqs. (28): which implies that the interaction kernel V (r, r ) between an isotropic dipole at r and one at r has the simple constant form for clamped boundaries, which is negative (representing an attractive interaction). The corresponding interaction for free boundaries is, which is positive (representing a repulsive interaction). Both the clamped and free boundary conditions lead to pairwise interactions between dipoles that are independent of the distance of separation between dipoles, and are inversely related to the system volume. They come with opposite signs and differ in strength by a factor of 2 but have the same generic behavior as seen above. The inverse volume scaling is responsible for the relation in Eq. (5) in the main text, whereas the sign of V 0 in Eq. (4) in the main text is determined according to the mechanical boundary condition by Eq. (30) and (31) above. Physically this kernel V(r,r') represents the interaction between a dipole and the images of the other dipoles 6 , since there are no direct elastic interactions between isotropic dipoles in an infinite medium. For a more general, anisotropic distribution of dipoles, the elastic interactions of dipoles in a finite 3D sphere can be expanded in the basis of spherical harmonics as 7 , where M l is negative (positive) if the boundaries are clamped (free), and scales inversely as the system volume.
The magnitude of M l decreases monotonically with l and we see from the form of Eq. (32) that the zeroth order contribution is dominant especially for dipoles far from the boundaries of the sphere. This justifies our consideration of an isotropic distribution of dipoles in this paper while also showing that higher order interactions are long-ranged and scale similarly with system volume; they are thus expected to give the same scaling relations for the mechanogen gradient vs. system size. Finally, we remark that since the mechanogen tensor is a scalar, it couples directly only to the trace of the strain or dipole tensors.  4) in 3D with only the gradient and quadratic terms (which are the dominant ones when the concentration is relatively large), in practice we see an intermediate region which better fits a pre-factor of unity, namely φ2 = a 2 /r 2 . The difference of a factor of two in the numerical pre-factor is probably a "crossover effect" arising from having to match the outer exponential solution with the inner power law behavior over numerical values of r that do not give a good enough separation of length scales. We have verified numerically that the pre-factor of the 1/r 2 term tends to increase from 1 to 2 as the separation of the scales r0 and λ(R) is increased, thereby providing further evidence that these are crossover effects. (Color online) Self-consistent solutions in 3D vs. radial coordinate rescaled with system size, R, showing approximate scaling:linear plots We solve the diffusion-degradation equation in 3D self-consistently, numerically for three different system sizes R1 = 20, 000a (blue), R2 = 10, 000a (red) and R3 = 5000a (green)and find the corresponding decay lengths λ(R)= 1000a, 520a and 280a respectively from the self-consistency condition λ(R) ∼ φ stated in the main text. This is a linear plot of the same solutions plotted on a semilog scale in Fig. (2) of the main text and exhibits approximate scaling as reported there.