# Mean-field interactions between living cells in linear and nonlinear elastic
matrices

Chaviva Sirote  Department of Biomedical Engineering, Tel Aviv University, Tel
Aviv 69978, Israel  Yair Shokef<EMAIL_ADDRESS>https://shokef.tau.ac.il
School of Mechanical Engineering, Tel Aviv University, Tel Aviv 69978, Israel
Sackler Center for Computational Molecular and Materials Science, Tel Aviv
University, Tel Aviv 69978, Israel  Center for Physics and Chemistry of Living
Systems, Tel Aviv University, Tel Aviv 69978, Israel

###### Abstract

Living cells respond to mechanical changes in the matrix surrounding them by
applying contractile forces that are in turn transmitted to distant cells. We
calculate the mechanical work that each cell performs in order to deform the
matrix, and study how that energy changes when a contracting cell is
surrounded by other cells with similar properties and behavior. We consider
simple effective geometries for the spatial arrangement of cells, with
spherical and with cylindrical symmetries, and model the presence of
neighboring cells by imposing zero-displacement at some distance from the
cell, which represents the surface of symmetry between neighboring cells. In
linear elastic matrices, we analytically study the dependence of the resulting
interaction energy on the geometry and on the stiffness and regulatory
behavior of the cells. For cells that regulate the active stress that they
apply, in spherical geometry, the deformation inside the cell is pure
compression thus the interaction depends only on their bulk modulus, while in
cylindrical geometries the deformation includes also shear and the interaction
depends also on their shear modulus. In nonlinear, strain stiffening matrices,
our numerical solutions and analytical approximations show how in the presence
of other cells, cell contraction is limited due to the divergence of the shear
stress.

##  I Introduction

Living cells embedded in, or adhered on an elastic environment transmit
mechanical forces through deformations of their surrounding solid medium
SchwarzRMP2013 . Thus their active contraction may be felt by non-contacting
cells through the propagation of stress and strain fields via this matrix.
Such mechanical interactions between cells and their surrounding environments,
as well as matrix-mediated cell-cell interactions are important for many
biological processes, such as stem-cell differentiation Engler2006 , wound
healing Poujade2007 , embryonic development Montrell2003 , cell division
Lesman2014 ; Abuhattum2015 , cancer metastasis Freidl2009 ; Gal2012 , and
cell-cell biochemical communication Gomez2019 ; Gomez2020 ; Jung2020 . Recent
experiments in synthetic setups or geometries enable to isolate and study
these mechanical interactions ReinhartKing2008 ; Winer2009 ; Buxboim2010 ;
Mohammadi2014 ; Lesman2015Interface ; Lesman2015Integrative .

Theoretically, cell contraction is customarily modeled by assuming that cells
generate active force dipoles Schwarz2002 , namely that they apply on the
elastic medium pairs of equal-magnitude and oppositely-pointing forces. The
interaction between distant such force dipoles may be described in terms of
the excess elastic energy stored in the medium in a situation with interacting
cells compared to the total energy of these cells without interactions.
Interaction energies between active force dipoles have been studied in the
past in the context of atoms adsorbed to surfaces Lau1977 , and more recently
in the cell mechanics context, with studies ranging from simple, linear
dipoles comprised of pairs of point forces Bischofs2003 ; Bischofs2004 ;
Bischofs2005 ; Bischofs2006  to more complex continuous shapes Golkov2019 .

To simplify the description of actively contractile cells, cells within a
three-dimensional medium Lesman2014 ; Lesman2015Integrative  have been modeled
as spherical force dipoles Shokef2012 ; Xu2015 ; BenYaakov2015 ; Wang2020 .
Namely, the mechanical activity of each cell is described by an isotropic
distribution of radial external forces that are applied on the surface of a
sphere within a three-dimensional elastic medium. This may be thought of as a
continuous collection of linear force dipoles distributed isotropically on the
surface of this sphere, with each such linear force dipole applying equal and
opposite forces at two opposing points on the sphere. This spherical geometry
implies spherical symmetry for the elastic fields around a cell, and thus
simplifies the partial differential equations of elastic equilibrium that
should be solved to a one-dimensional or ordinary differential equation for
the radial displacement as a function of the radial coordinate.

Along these lines, interactions between biological cells have been modeled by
imposing a zero-displacement boundary condition at some distance from this
spherical cell BenYaakov2015 . The justification for this is to first assume
that the spatial distribution of cells in their elastic environment is on a
periodic array. In such case, by symmetry, the normal displacement on the
surface of each periodic unit-cell of this array should vanish. A reasonable
mean-field approximation would then be to replace the actual polyhedral
periodic unit-cell with a spherical unit-cell with zero displacement on its
boundary, which results in a simple one-dimensional geometry also for the
study of cell-cell interactions. Note that the periodic unit-cell of an array
of biological cells and the spherical domain approximating this unit-cell
should not be confused with the biological cells.

Cell-cell elastic interactions are similar, but not identical to electrostatic
interactions between induced electric dipoles Israelachvili1992 ; Schwarz2002
; Bischofs2004 ; SchwarzRMP2013 ; Golkov2017 , and the interaction energy
scales algebraically with distance. Interestingly, the sign of the interaction
energy is set by the homeostatic behavior of the cells; cells regulating the
displacement on their boundary repel each other, and cells regulating the
force that they apply are mechanically attracted to each other BenYaakov2015 .
Elaborate solutions, which include the interference between the angular
dependence of the deformation fields around two spheres recover these mean-
field results, both in terms of the sign and of the exponents of the power-law
decay in the magnitude of the interaction energy Golkov2017 .

In this Paper we use the aforementioned mean-field approach to study
additional geometries, namely not only spherical symmetry, but also two simple
cylindrical setups Sopher2018  that relate to cells on substrates or to
elongated cells in a three-dimensional matrix. In doing that, we provide
further insight on matrix-mediated elastic interactions in more complex
geometries, for which analytical results are not known. Furthermore, we
include in our calculations the stiffness of the interior the cell, and not
only that of the surrounding matrix. Finally, and most importantly, we extend
our analysis to cells in nonlinear, strain-stiffening media Shokef2012 .

##  II Model and Interaction Energy

We analyze two geometries, one in which the contracting cell is assumed to be
a sphere of radius $R_{0}$, and the other in which it is assumed to be a
cylinder of radius $R_{0}$. In both geometries, we model the presence of
neighboring cells by imposing a rigid wall with zero normal displacement at
some larger radius $R_{1}$, see Fig. 1. This radius represents half of the
distance to nearby cells, since in realistic geometries, by symmetry the
normal displacement should vanish at that position. In our model, the
biological cell has linear elastic response characterized by shear modulus
$\mu_{c}$ and bulk modulus $K_{c}$. In Sec. III we will assume that the matrix
is also linearly elastic with shear modulus $\mu_{m}$ and bulk modulus
$K_{m}$, while in Sec. IV we will include nonlinear strain-stiffening of the
matrix.

Figure 1: Schematic of a single cell in an infinite matrix (a) and in a
bounded region (b). The cell may be a sphere (c), a flat cylinder (d) for
plane stress, or a long cylinder (e) for plane strain.

The elastic energy stored in the cell and in the matrix together is equal to
the mechanical work performed by the active, or external forces, which act
only at the surface of the cell. Here, these are radial forces that act at the
radius $R_{0}$ and generate a radial displacement $u_{0}$ there. We assume
that in the absence of these active forces, the elastic medium is relaxed, and
hence we calculate the work done on the matrix by integrating over the
adiabatic process of building up the displacement $u_{0}$, starting from this
relaxed state Golkov2019

$E=\int_{0}^{u_{0}}g\tau(w)dw.$ (1)

Here $w$ denotes the displacement at $R_{0}$ during this process, and
$\tau(w)=\sigma_{c}(R_{0})-\sigma_{m}(R_{0})$ is the active stress applied by
the cell, namely the discontinuity in the radial component of the stress
tensor on the surface of the cell when the displacement there is equal to $w$.
We use $\sigma_{c}$ and $\sigma_{m}$ to denote the radial component
$\sigma_{rr}$ of the stress tensor in the cell and in the matrix,
respectively. The geometrical factor relating stress to force is equal to
$g=4\pi R_{0}^{2}$ in the spherical geometry. In the cylindrical geometry we
use $g=2\pi R_{0}$ to obtain the force per unit length of the cylinder, and
subsequently the energy per unit length. Since much of our analysis and
discussion below is performed simultaneously for the two geometries, for
brevity we will use the term energy and the symbol E, for energy in the
spherical geometry and for energy per unit length in the cylindrical geometry.
If the medium is linearly elastic (Sec. III), there is a linear
proportionality between displacement and active stress, $\tau(w)\propto w$,
and Eq. (1) reduces to $E=\frac{1}{2}g\sigma_{0}u_{0}$, where
$\sigma_{0}=\tau(u_{0})$ is the active stress, or stress difference on the
surface of the cell when the displacement there is equal to $u_{0}$
BenYaakov2015 ; Golkov2017 . However, for nonlinear media (Sec. IV), the
entire dependence of $\tau$ on $w$ up to the actual displacement $u_{0}$ is
required in order to perform the integration in Eq. (1).

We define the interaction energy as the added work that the cell has to
perform due to the presence of neighboring cells. Within our mean-field
approximation, this is equal to the difference

$\Delta E=E_{1}-E_{\infty}$ (2)

between the elastic energy $E_{1}$ stored in the medium when there is a zero-
displacement boundary condition at $R_{1}$ (Fig. 1b) and the elastic energy
$E_{\infty}$ of the same cell in an unbounded matrix (Fig. 1a). Care should be
given to the meaning of this placement of the same cell in two different
mechanical environments. Obviously when identical cells have different
environments, they behave differently Engler2006 ; Brown1998 ; Freyman2002 ;
Saez2005 ; Ghibaudo2008 ; De2008 ; He2014 . In our case, the confining
geometry alters the relation between the active stress $\tau$ and the
displacement $u_{0}$. Thus cells in the confined and in the unbounded
geometries cannot simultaneously have both the same active stress and the same
deformation. We will consider two extreme scenarios: i) displacement
regulation, for which cells are assumed to be biologically programmed to
generate a given displacement $u_{0}$, regardless of what active stress they
need to apply in order to reach that deformation, and ii) stress regulation,
for which cells apply a given active stress $\sigma_{0}$, regardless of the
displacement they manage to generate.

For cells regulating the displacement on their surface, the far-field zero-
displacement boundary condition imposed by the presence of neighboring cells
causes cells to apply a higher active stress in order to generate a given
contraction. Thus we expect the interaction energy to be positive,
representing repulsion between cells. For stress regulation, on the other
hand, a smaller displacement will be generated in the presence of other cells,
thus leading to a negative interaction energy, or to attraction BenYaakov2015
; Golkov2017 . We consider cells that are embedded in a solid surrounding,
thus they cannot easily move due to these repulsive or attractive
interactions. However, we suggest that these interactions are the mechanical
cause for instance for cells to send protrusions one toward each other or
rather away from each other Mohammadi2014  or for cells to generate bands of
densification and alignment of the network of fibers around them Gomez2019 ;
Natan2021 , and also for migratory behavior of cells on substrates
ReinhartKing2008 . Alternatively, our results could be used to explain when do
cells regulate their displacement and when do they regulate their active
stress, depending on the energetic cost for doing so.

Our strategy for calculating the interaction energy for the different cases
described above is as follows: We first assume a displacement $u_{0}$ on the
boundary of the cell, and solve for the displacement field. Using the
resulting strain, we calculate the stress on the boundary of the cell, and
from that the active stress $\sigma_{0}$ that the cell has to apply to
generate the displacement $u_{0}$. We then use these to calculate the elastic
energy stored in the deformation. For displacement regulation we subtract the
resulting energies in the bounded and in the unbounded geometries to get the
interaction energy, Eq. (2) at a given displacement $u_{0}$. For stress
regulation, we combine our results for $E$ vs. $u_{0}$ and for $\sigma_{0}$
vs. $u_{0}$ to obtain $E$ vs. $\sigma_{0}$, and then use Eq. (2) to get the
interaction energy at a given active stress $\sigma_{0}$. This procedure has
been applied previously BenYaakov2015  for an empty spherical cell in a
linearly elastic medium, with the focus on the asymptotic behavior when cells
are far apart from each other, which translates to $R_{1}\gg R_{0}$ in our
model. Here, we provide the full dependence on distance, and also: i) include
a passive, elastic rigidity of the interior of the cell, ii) solve for
additional cases with cylindrical symmetry, and iii) study the effects of
nonlinear elastic response of the matrix.

##  III Linear Medium

Our spherical and cylindrical geometries both assume radial symmetry, such
that the displacement vector $\vec{u}(\vec{r})=u(r)\hat{r}$ is given only by
the radial component $u(r)$, which depends only on the radial coordinate $r$.
In this section we assume that the mechanical response of both the cell and
the surrounding matrix is linearly elastic. This leads to linear ordinary
differential equations for $u(r)$, which we analytically solve. Thus we obtain
exact analytical expressions for the interaction energy in the different
geometries. More complicated nonlinear material models may be considered for
the cell and for the matrix, for instance the nonlinear elastic response of
the matrix as described in Sec. IV below. Such models result in nonlinear
differential equations that need to be solved in order to obtain the
displacement field. However, the spherical and cylindrical symmetries of our
geometrical setups reduces the general three-dimensional partial differential
equations of nonlinear elasticity to one-dimensional ordinary differential
equations, which, if needed, may be solved numerically much more easily.

###  III.1 Spherical Cell

For spherical symmetry (Fig. 1c), mechanical equilibrium implies Bower2009

$\frac{d^{2}u}{dr^{2}}+\frac{2}{r}\frac{du}{dr}-\frac{2u}{r^{2}}=0.$ (3)

The general solution to this equation is:

$u(r)=Cr+\frac{D}{r^{2}},$ (4)

where the constants $C$ and $D$ are set by the boundary conditions. For
spherical symmetry, the radial stress is given by Bower2009

$\sigma_{rr}=\frac{4}{3}\mu\left(\frac{du}{dr}-\frac{u}{r}\right)+K\left(\frac{du}{dr}+\frac{2u}{r}\right),$
(5)

which for the general solution (4) may be written as

$\sigma_{rr}=3KC-\frac{4\mu D}{r^{3}}.$ (6)

Thus for spherical symmetry, the $u\propto r$ term in (4) represents pure
compression, while the $u\propto 1/r^{2}$ term represents pure shear.

Inside the cell ($0<r<R_{0}$) the boundary conditions are $u(0)=0$ and
$u(R_{0})=u_{0}$, thus the displacement is given by

$u_{c}=\frac{u_{0}r}{R_{0}},$ (7)

which contains only compression, and the stress in the cell is uniform and
equal to

$\sigma_{c}=3K_{c}\frac{u_{0}}{R_{0}},$ (8)

which depends only on the bulk modulus of the cell and not on its shear
modulus.

In the surrounding matrix ($r>R_{0}$), the displacement at the cell border is
$u(R_{0})=u_{0}$. For the unbounded case, the second boundary condition is
that $u=0$ at $r=\infty$, leading to

$u_{m,\infty}=\frac{u_{0}R_{0}^{2}}{r^{2}},$ (9)

which contains only shear. The stress in the matrix in this case is thus

$\sigma_{m,\infty}=-\frac{4\mu_{m}u_{0}R_{0}^{2}}{r^{3}}=-4\mu_{m}\frac{u_{0}}{R_{0}},$
(10)

where the second equality is obtained by setting the position to be on the
surface of the cell, $r=R_{0}$. The total work that an isolated cell performs
in order to deform both itself and the matrix that surrounds it is thus

$E_{\infty}=2\pi u_{0}^{2}R_{0}(3K_{c}+4\mu_{m}).$ (11)

For the bounded case, the second boundary condition is now $u(R_{1})=0$, and
the displacement in the matrix is

$u_{m,1}=\frac{u_{0}R_{0}^{2}}{R_{1}^{3}-R_{0}^{3}}\left(\frac{R_{1}^{3}}{r^{2}}-r\right),$
(12)

which contains both shear and compression. In this case, the stress in the
matrix is

$\sigma_{m,1}=-\frac{4\mu_{m}R_{1}^{3}+3K_{m}R_{0}^{3}}{R_{1}^{3}-R_{0}^{3}}\frac{u_{0}}{R_{0}},$
(13)

and the energy is

$E_{1}=2\pi
u_{0}^{2}R_{0}\left[3K_{c}+\frac{4\mu_{m}R_{1}^{3}+3K_{m}R_{0}^{3}}{R_{1}^{3}-R_{0}^{3}}\right].$
(14)

####  III.1.1 Displacement Regulation

If the cell generates the same displacement $u_{0}$ on its boundary regardless
of its mechanical environment, the interaction energy is just the difference
between the energy (14) when it is bounded and the energy (11) when it is
unbounded, namely

$\Delta E=E_{1}-E_{\infty}=2\pi
u_{0}^{2}R_{0}^{4}\frac{4\mu_{m}+3K_{m}}{R_{1}^{3}-R_{0}^{3}}.$ (15)

This result depends on both the bulk and shear moduli of the matrix because
both compression and shear are generated in the matrix. It does not, however,
depend on the elastic moduli of the cell since for displacement regulation,
the cell deforms in the exact same way in the bounded and in the unbounded
geometries, thus the same energy is stored in the cell for both geometries. As
predicted above, for displacement regulation, the interaction energy is
positive, representing repulsion between cell. For large separation between
cell ($R_{1}\gg R_{0}$), the interaction energy decays as $1/R_{1}^{3}$, which
is consistent with a mapping to interactions between induced electric dipoles
BenYaakov2015 ; Golkov2017 .

####  III.1.2 Stress Regulation

If the cell generates the same active stress $\sigma_{0}$ both when it is
bounded and when it is unbounded, we need to first express the energies
$E_{1}$ and $E_{\infty}$ in terms of $\sigma_{0}$ rather than $u_{0}$, and
only then subtract them to get the interaction energy. Using Eqs. (10) and
(8), for an unbounded cell, the displacement on its boundary is given by

$u_{0,\infty}=\frac{\sigma_{0}R_{0}}{3K_{c}+4\mu_{m}},$ (16)

And from Eq. (11), the energy in this case is given by

$E_{\infty}=\frac{2\pi\sigma_{0}^{2}R_{0}^{3}}{3K_{c}+4\mu_{m}}.$ (17)

Similarly, from Eqs. (13) and (8), the displacement on the boundary of a
bounded cell is

$u_{0,1}=\frac{\sigma_{0}R_{0}(R_{1}^{3}-R_{0}^{3})}{4\mu_{m}R_{1}^{3}+3K_{m}R_{0}^{3}+3K_{c}\left(R_{1}^{3}-R_{0}^{3}\right)},$
(18)

and the energy (14) in this case is

$E_{1}=\frac{2\pi\sigma_{0}^{2}R_{0}^{3}(R_{1}^{3}-R_{0}^{3})}{4\mu_{m}R_{1}^{3}+3K_{m}R_{0}^{3}+3K_{c}\left(R_{1}^{3}-R_{0}^{3}\right)}.$
(19)

By subtracting (17) from (19) we get that for stress regulation the
interaction energy is

$\Delta
E=-\frac{2\pi\sigma_{0}^{2}R_{0}^{6}(4\mu_{m}+3K_{m})}{\left(4\mu_{m}+3K_{c}\right)\left[4\mu_{m}R_{1}^{3}+3K_{m}R_{0}^{3}+3K_{c}\left(R_{1}^{3}-R_{0}^{3}\right)\right]}.$
(20)

As discussed above, for stress regulation, the interaction energy is negative,
representing attraction between cells. Moreover, here the contraction of the
cell differs between the bounded and unbounded geometries, therefore the
interaction energy depends on the stiffness of the cell. However, for
spherical symmetry the deformation of the cell includes only compression, thus
the interaction energy depends on the bulk modulus of the cell but not on its
shear modulus. The $1/R_{1}^{3}$ decay in the magnitude of the interaction
energy for large distances between cells is similar to that found above for
displacement regulation.

###  III.2 Cylindrical Cell

For cylindrical symmetry, the boundary conditions at the ends of the cylinder
affect the stress inside it. If the cylinder is prevented from deforming in
the axial direction, the axial strain vanishes $\epsilon_{zz}=0$, this is
termed plane strain, and the radial stress reads Bower2009

$\sigma_{rr}^{A}=\frac{2}{3}\mu\left(2\frac{du}{dr}-\frac{u}{r}\right)+K\left(\frac{du}{dr}+\frac{u}{r}\right).$
(21)

In order to maintain this state, external axial forces should be applied at
the ends of the cylinder. Note that there is no displacement there, hence
these forces do not perform any work, and thus do not contribute to the
energy. If such forces are not applied, the axial stress is zero
$\sigma_{zz}=0$, this is termed plane stress, and the radial stress is given
by Bower2009

$\sigma_{rr}^{E}=\frac{2}{3}\mu\left(2\frac{du}{dr}-\frac{u}{r}\right)+\frac{1-4\nu^{2}}{1-\nu^{2}}K\left(\frac{du}{dr}+\frac{u}{r}\right),$
(22)

with the Poisson ratio given by

$\nu=\frac{3K-2\mu}{2(3K+\mu)}.$ (23)

Also for plane stress no work is done due to the axial deformation since no
forces are applied in that direction. The simple geometries considered here
are not meant to exactly describe the shapes of specific biological cells and
their environments. However, we suggest that plane stress (Fig. 1d) may be
more relevant for thin cells in a monolayer or for cells on a substrate, while
plane strain (Fig. 1e) may be more relevant for elongated cells such as
neurons, in a three dimensional environment.

Regardless of the aforementioned boundary conditions in the axial direction,
cylindrical symmetry implies that the radial displacement satisfies the
following differential equation, which ensures mechanical equilibrium
Bower2009

$\frac{d^{2}u}{dr^{2}}+\frac{1}{r}\frac{du}{dr}-\frac{u}{r^{2}}=0.$ (24)

The general solution to this equation reads:

$u(r)=Cr+\frac{D}{r},$ (25)

where the constants $C$ and $D$ are set by the boundary conditions in the
radial direction. We now generalize the two cases for the axial boundary
conditions, plane strain (21) and plane stress (22) to write the radial stress
for both cases as

$\sigma_{rr}=\frac{2}{3}\mu\left(2\frac{du}{dr}-\frac{u}{r}\right)+\hat{K}\left(\frac{du}{dr}+\frac{u}{r}\right),$
(26)

with an effective bulk modulus given by

$\hat{K}=\left\\{\begin{array}[]{cc}K&\rm{plane\>strain}\\\
\frac{1-4\nu^{2}}{1-\nu^{2}}K&\rm{plane\>stress}\end{array}\right.$ (27)

Substituting the general solution for the displacement (25) in (26) enables to
generally write the stress as

$\sigma_{rr}=2\mu\left(\frac{C}{3}-\frac{D}{r^{2}}\right)+2\hat{K}C.$ (28)

From this we see that the solution $u=D/r$ contains only shear. However, as
opposed to the spherical case discussed above, in cylindrical geometry, the
solution $u=Cr$ contains both shear and compression.

Similarly to the procedure employed above for a spherical cell, we will now
solve for the displacement field $u(r)$, from that we will obtain the relation
between the displacement $u_{0}$ on the cell boundary and the active stress
$\sigma_{0}$ that the cell applies, and from that we will calculate the
elastic energy for a bounded and for an unbounded cell, and subsequently the
interaction energy for displacement regulation and for stress regulation.

Inside the cell, the displacement has the same form as for the spherical
geometry,

$u_{c}=\frac{u_{0}r}{R_{0}}.$ (29)

Here, by Eq. (28), the stress inside the cell is given by

$\sigma_{c}=\left(\frac{2}{3}\mu_{c}+2\hat{K}_{c}\right)\frac{u_{0}}{R_{0}}.$
(30)

Outside the cell, for the unbounded cylindrical geometry, the displacement is

$u_{m,\infty}=\frac{u_{0}R_{0}}{r},$ (31)

which leads to the following radial stress on the surface of the cell
($r=R_{0}$)

$\sigma_{m,\infty}=-2\mu_{m}\frac{u_{0}}{R_{0}},$ (32)

and subsequently to the energy per unit length

$E_{\infty}=2\pi u_{0}^{2}\left(\frac{\mu_{c}}{3}+\hat{K}_{c}+\mu_{m}\right).$
(33)

It is interesting to note that $E_{\infty}$ does not depend on the cell
radius, $R_{0}$. For the bounded geometry, the displacement is given by

$u_{m,1}=\frac{u_{0}R_{0}}{R_{1}^{2}-R_{0}^{2}}\left(\frac{R_{1}^{2}}{r}-r\right),$
(34)

thus the stress on the surface of the cell reads

$\sigma_{m,1}=-\frac{2u_{0}R_{0}}{R_{1}^{2}-R_{0}^{2}}\left[\mu_{m}\left(\frac{1}{3}+\frac{R_{1}^{2}}{R_{0}^{2}}\right)+\hat{K}_{m}\right],$
(35)

and the energy per unit length is

$\displaystyle E_{1}=$ $\displaystyle\pi
u_{0}^{2}\left\\{\frac{2}{3}\mu_{c}+2\hat{K}_{c}+\right.$
$\displaystyle\left.\left[\frac{2}{3}\mu_{m}\left(1+3\frac{R_{1}^{2}}{R_{0}^{2}}\right)+2\hat{K}_{m}\right]\frac{R_{0}^{2}}{R_{1}^{2}-R_{0}^{2}}\right\\}.$
(36)

As for the spherical case, for fixed displacement, we obtain the interaction
energy by subtracting (33) from (36), leading to

$\Delta E=\frac{2}{3}\pi
u_{0}^{2}R_{0}^{2}\frac{4\mu_{m}+3\hat{K}_{m}}{R_{1}^{2}-R_{0}^{2}}.$ (37)

Similarly to the procedure described above for spherical symmetry, also in
cylindrical symmetry for regulation of the active stress, we use (30), (32)
and (35) to obtain $u_{0}$ in terms of $\sigma_{0}$ for the bounded and for
the unbounded cases, and from that obtain the interaction energy,

$\Delta
E=-\frac{\pi}{6}\sigma_{0}^{2}R_{0}^{4}\frac{\left(4\mu_{m}+3\hat{K}_{m}\right)}{\left(\mu_{m}+\frac{\mu_{c}}{3}+\hat{K}_{c}\right)\left[\left(\frac{\mu_{m}}{3}+\hat{K}_{m}\right)R_{0}^{2}+\mu_{m}R_{1}^{2}+\left(\hat{K}_{c}+\frac{\mu_{c}}{3}\right)\left(R_{1}^{2}-R_{0}^{2}\right)\right]}.$
(38)

As for the spherical geometry, also for the cylindrical geometry we get for
displacement regulation a repulsive interaction (37) which depends only on the
elastic moduli of the matrix, while for stress regulation we obtain attractive
interaction (38), which depends also on the elastic moduli of the cell. Here
there is both shear and compression in the cell, thus (38) depends on both
$\mu_{c}$ and $K_{c}$. For large separations between cells, $R_{1}\gg R_{0}$,
the interaction energy in the cylindrical geometry decays as $1/R_{1}^{2}$, as
opposed to $1/R_{1}^{3}$ in spherical geometry, which is consistent with the
dimensionality reduction between the two cases.

##  IV Nonlinear Medium

###  IV.1 Nonlinear Material Model

The analysis presented in Sec. III assumes linear elasticity and is hence
valid only for extremely small deformations. The extracellular matrix is a gel
of crosslinked polymers, which exhibits nonlinear mechanical response; its
response to shear is linear at small stress, while for higher stress, the
differential shear modulus crosses over to a power law increase with shear
stress, $G(\sigma)\propto\sigma^{3/2}$ Gardel2004 ; Storm2005 ; Vader2009 ;
Broedersz2014 . Due to their near-incompressibility, biopolymer gels are
usually compressed much less than they are sheared, thus we shall include
nonlinearity only in their response to shear. Specifically, we describe the
nonlinear medium surrounding our cells using the following energy density
function Shokef2012

$W=\frac{\mu}{2}\left\\{\left[1-b\left(\bar{I}_{1}-3\right)\right]^{-1}-1\right\\}+\frac{K}{2}\left(J-1\right)^{2},$
(39)

where $\mu$ and $K$ are the shear and bulk moduli in the linear regime, and
$b$ is a dimensionless parameter quantifying the strength of the nonlinearity.
The second term describes compression in a neo-Hookean manner with
$J=\det(\mathbf{F})$ quantifying the compression, and $F_{i,j}=\partial
x_{i}/\partial X_{j}$ the deformation gradient tensor, where $\vec{X}$ is the
reference position and $\vec{x}$ the deformed position. The first term in (39)
describes shear with $\bar{I}_{1}=\mathrm{tr}(\mathbf{B})/J^{2/3}$ the
normalized first eigenvalue of the left Cauchy-Green strain tensor,
$\mathbf{B}=\mathbf{F}\mathbf{F}^{T}$. Most generally, the Cauchy stress
tensor is derived from the energy density as

$\sigma_{ij}=\frac{1}{J}F_{ik}\frac{\partial W}{\partial F_{kj}}.$ (40)

However, for an energy density function which depends only on $J$ and $I_{1}$,
this reduces to Bower2009

$\boldsymbol{\sigma}=\frac{2}{J^{5/3}}\frac{\partial
W}{\partial\bar{I}_{1}}\left({\bf B}-\frac{I_{1}}{3}{\bf
1}\right)+\frac{\partial W}{\partial J}{\bf 1},$ (41)

with ${\bf 1}$ the unit tensor. For our model (39), this yields

$\boldsymbol{\sigma}=\frac{\mu}{J^{5/3}}\left[1-b\left(\bar{I}_{1}-3\right)\right]^{-2}\left({\bf
B}-\frac{I_{1}}{3}{\bf 1}\right)+K(J-1){\bf 1}.$ (42)

Figure 2: Stress stiffenning for simple shear: a) Stress vs strain is linear
for small strain and then diverges at $\gamma_{a}=1/\sqrt{b}$. b) differential
shear modulus is constant for small stress and then increases as
$G\propto\sigma^{3/2}$ for stress larger than $\mu\gamma_{a}$.

For simple shear, the deformation may be written as $x=X+\gamma Z$, $y=Y$,
$z=Z$, thus Eq. (42) gives a shear stress of
$\sigma_{xz}=\mu\gamma\left(1-b\gamma^{2}\right)^{-2}$, resulting in the
differential shear modulus $G\equiv
d\sigma/d\gamma=\mu\left(1+3b\gamma^{2}\right)\left(1-b\gamma^{2}\right)^{-3}$.
In the limit of small shear $\left(\gamma\ll 1\right)$, we obtain linear
response $\sigma\approx\mu\gamma$ with a constant differential shear modulus
equal to $\mu$. For large shear, the shear stress $\sigma_{xz}$ diverges as
the strain $\gamma$ approaches a maximal possible strain
$\gamma_{a}=1/\sqrt{b}$. In that limit of $\sigma\gg\mu\gamma_{a}$, the
differential shear modulus may be approximated as $G\approx
4\mu\left[\sigma/\left(\mu\gamma_{a}\right)\right]^{3/2}$, which exhibits the
power law stiffening $G\propto\sigma^{3/2}$ characteristic of biopolymer gels,
see Fig. 2.

Below, we will study the effects of this power-law stiffening on matrix
mediated elastic interactions between cells. We will consider the same
geometrical model studied above, of a symmetric contracting cell, surrounded
by an effective wall, which provides a mean-field description for the effect
of neighboring cells. For simplicity, we will restrict ourselves to a
spherical cell. However, the extensions described above for for cylindrical
cells are straightforward.

In order to simplify the analysis, we will restrict ourselves to the small
deformations limit, which is relevant for strongly nonlinear materials ($b\gg
1$) that have a very small maximal shear strain $\gamma_{a}\ll 1$. In such
cases, the deformations will always be small $\gamma\ll 1$, and we will have
strong nonlinear effects even at small deformations. In our spherical case,
this implies $u/r,du/dr\ll 1$. We will obviously keep terms linear in $u$.
However, since for simple shear $b=1/\gamma_{a}^{2}$, for consistency, we will
keep also terms proportional to $bu^{2}$, and will neglect only terms that are
higher order than that.

For spherical symmetry, the deformation is described by the difference
$u(R)=r(R)-R$ between the deformed radius $r$ and the reference radius $R$.
The compression and shear components of the deformation are generally given by
Bower2009

$\displaystyle J$ $\displaystyle=$
$\displaystyle\left(1+\frac{u}{R}\right)^{2}\left(1+\frac{du}{dR}\right),$
$\displaystyle\bar{I}_{1}$ $\displaystyle=$
$\displaystyle\left(1+\frac{u}{R}\right)^{-4/3}\left(1+\frac{du}{dR}\right)^{4/3}$
(43) $\displaystyle+$ $\displaystyle
2\left(1+\frac{u}{R}\right)^{2/3}\left(1+\frac{du}{dR}\right)^{-2/3}.$

For small displacements ($u\ll R$), these are approximated as

$\displaystyle J$ $\displaystyle\approx$ $\displaystyle
1+2\frac{u}{r}+\frac{du}{dr},$ $\displaystyle\bar{I}_{1}$ $\displaystyle-$
$\displaystyle 3\approx\frac{4}{3}\left(\frac{du}{dr}-\frac{u}{r}\right)^{2},$
(44)

where we have also used the small deformation approximation to replace $R$
with $r$.

Force balance $\nabla\cdot\boldsymbol{\sigma}=0$ implies

$\frac{d\sigma_{rr}}{dr}+\frac{2}{r}\left(\sigma_{rr}-\sigma_{\theta\theta}\right)=0.$
(45)

In spherical coordinates, the general expression for the stress tensor (41)
gives

$\displaystyle\sigma_{rr}=\frac{4}{3J^{5/3}}\frac{\partial
W}{\partial\bar{I}_{1}}\left[\left(\frac{du}{dr}\right)^{2}+2\frac{du}{dr}-\left(\frac{u}{r}\right)^{2}-2\frac{u}{r}\right]+\frac{\partial
W}{\partial J},$
$\displaystyle\sigma_{rr}-\sigma_{\theta\theta}=\frac{2}{J^{5/3}}\frac{\partial
W}{\partial\bar{I}_{1}}\left[\left(\frac{du}{dr}\right)^{2}+2\frac{du}{dr}-\left(\frac{u}{R}\right)^{2}-2\frac{u}{r}\right],$
(46)

where for our strain energy density (39),

$\frac{\partial
W}{\partial\bar{I}_{1}}=\frac{\mu}{2}\left[1-b(\bar{I}_{1}-3)\right]^{-2},$
(47)

and

$\frac{\partial W}{\partial J}=K\left(\frac{du}{dr}+\frac{2u}{r}\right).$ (48)

By substituting the partial differentials of the energy density function
(47,48) in the expression (46) for the stress tensor, we may write the
condition for mechanical equilibrium (45) as

$\displaystyle\left\\{K+\frac{4\mu}{3}\left[1-b\left(\bar{I}_{1}-3\right)\right]^{-2}\right\\}\left(\frac{d^{2}u}{dr^{2}}+\frac{2}{r}\frac{du}{dr}-\frac{2u}{r^{2}}\right)+\frac{8\mu}{3}b\left[1-b\left(\bar{I}_{1}-3\right)\right]^{-3}\frac{d\bar{I}_{1}}{dr}\left(\frac{du}{dr}-\frac{u}{r}\right)=0.$
(49)

After multiplying by $\left[1-b\left(\bar{I}_{1}-3\right)\right]^{3}$ we write

$\displaystyle\left[1-b\left(\bar{I}_{1}-3\right)\right]\left\\{K\left[1-b\left(\bar{I}_{1}-3\right)\right]^{2}+\frac{4\mu}{3}\right\\}\left(\frac{d^{2}u}{dr^{2}}+\frac{2}{r}\frac{du}{dr}-\frac{2u}{r^{2}}\right)+\frac{8\mu}{3}b\frac{d\bar{I}_{1}}{dr}\left(\frac{du}{dr}-\frac{u}{r}\right)=0.$
(50)

To bring (50) to dimensionless form, we normalize the radial coordinate by the
cell radius $\tilde{r}=r/R_{0}$ and the displacement by the displacement on
the cell boundary $\tilde{u}=u/u_{0}$. We quantify the nonlinearity by
$A=b\left(u_{0}/R_{0}\right)^{2}=\left[\left(u_{0}/R_{0}\right)/\gamma_{a}\right]^{2}$,
which is a measure of the ratio between the typical scale $u_{0}/R_{0}$ of the
strain on the cell surface and the maximal shear strain
$\gamma_{a}=1/\sqrt{b}$ that the material can sustain. Finally, we express the
ratio of bulk modulus to linear shear modulus in terms of the Poisson ratio
(23), and eventually write the equilibrium condition as:

$\displaystyle\left(1-A\tilde{I}_{3}\right)\left[\frac{1+\nu}{2(1-2\nu)}\left(1-A\tilde{I}_{3}\right)^{2}+1\right]\left(\frac{d^{2}\tilde{u}}{d\tilde{r}^{2}}+\frac{2}{\tilde{r}}\frac{d\tilde{u}}{d\tilde{r}}-2\frac{\tilde{u}}{\tilde{r}^{2}}\right)+2A\frac{d\tilde{I_{3}}}{d\tilde{r}}\left(\frac{d\tilde{u}}{d\tilde{r}}-\frac{\tilde{u}}{\tilde{r}}\right)=0,$
(51)

where

$\displaystyle\tilde{I}_{3}=\left(\frac{R_{0}}{u_{0}}\right)^{2}\left(\bar{I}_{1}-3\right)=\frac{4}{3}\left(\frac{d\tilde{u}}{d\tilde{r}}-\frac{\tilde{u}}{\tilde{r}}\right)^{2}$
(52)

is a dimensionless measure of the shear in the matrix.

Equation (51) is a nonlinear second order ordinary differential equation for
the dimensionless displacement as a function of the dimensionless radius,
$\tilde{u}\left(\tilde{r}\right)$. We obtain $\tilde{u}\left(\tilde{r}\right)$
by numerically solving the boundary value problem using our two boundary
conditions: given displacement on the cell border $\tilde{u}(1)=1$, and zero
displacement, either at infinity for the unbounded geometry
$\tilde{u}(\infty)=0$ or at $\tilde{R}_{1}=R_{1}/R_{0}$ for the bounded
geometry $\tilde{u}\left(\tilde{R}_{1}\right)=0$.

Clearly, in the limit $A=0$ of vanishing nonlinearity, Eq. (51) reduces to Eq.
(3), which describes a linearly elastic medium. In this linear limit, there
are two solutions, $\tilde{u}=\tilde{r}$ and $\tilde{u}=1/\tilde{r}^{2}$, and
as shown in Sec. III, a linear combination of the two can match any given set
of boundary conditions. More interestingly, in this spherical geometry, the
solution $\tilde{u}=\tilde{r}$ that contains only compression without shear,
leads to $\tilde{I}_{3}=0$ and thus solves Eq. (51) for arbitrary strength $A$
of the nonlinearity. This solution is consistent with the first boundary
condition at the cell border $\tilde{u}(1)=1$, but it increases in magnitude
with $\tilde{r}$ and thus cannot satisfy the second, zero displacement
boundary condition at $\tilde{r}=\infty$ or at $\tilde{r}=\tilde{R}_{1}$.

###  IV.2 Unbounded Cell

For an unbounded cell, the dimensionless displacement
$\tilde{u}\left(\tilde{r}\right)$ depends only on the dimensionless strength
$A$ of the nonlinearity and on the Poisson ratio $\nu$. The behavior of an
unbounded cell in this nonlinear material model depends very weakly on the
Poisson ratio Shokef2012 , thus we will not consider this dependence here, and
all the numerical results we present below are for $\nu=0.4$. As shown in Fig.
3a, for $A=0$ the response is completely linear, and the displacement decays
according to the prediction of linear elasticity,
$\tilde{u}\left(\tilde{r}\right)=1/\tilde{r}^{2}$, or
$u(r)=u_{0}\left(R_{0}/r\right)^{2}$. As $A$ increases, an increasingly larger
region near the cell responds nonlinearly, approaching the shear-less solution
$\tilde{u}\left(\tilde{r}\right)=\tilde{r}$, or $u(r)=u_{0}r/R_{0}$, while far
enough from the cell $\tilde{r}\gg 1$, deformations decay and the matrix
restores its linear response, leading to a $\tilde{u}\propto 1/\tilde{r}^{2}$
decay of displacements in the far field. The far-field behavior may be written
as $u(r)=u_{\rm eff}\left(R_{0}/r\right)^{2}$, or
$\tilde{u}\left(\tilde{r}\right)=\tilde{u}_{\rm eff}/\tilde{r}^{2}$, where
$\tilde{u}_{\rm eff}=u_{\rm eff}/u_{0}$ describes the amplification of the
effective displacement on the cell border as felt at long distances. We will
denote by $\tilde{R}_{*}\equiv R_{*}/R_{0}$ the position of the maximum of
$\tilde{u}(\tilde{r})$. This position roughly divides space into the nonlinear
near-field compression-dominated regime at $\tilde{r}<\tilde{R}_{*}$ and the
linearly-elastic far-field shear-dominated regime at
$\tilde{r}>\tilde{R}_{*}$.

Figure 3: Normalized displacement $\tilde{u}=u/u_{0}$ vs. normalized radius
$\tilde{r}=r/R_{0}$ in spherical nonlinear medium for the (a) unbounded
geometry and (b) geometry bounded at
$\tilde{R}_{1}=R_{1}/R_{0}=10,100,250,500,1000$. Colors mark the dimensionless
strength $A=bu_{0}^{2}/R_{0}^{2}$ of the nonlinearity, as indicated in the
legend. All results are for Poisson ratio $\nu=0.4$. In both geometries, for
increasing nonlinearity, the near-field behavior tends toward the shear-less
solution $\tilde{u}=\tilde{r}$ represented by the blue dotted line, and may be
approximated by Eq. (53), which is drawn in a black dashed line only for the
maximal $A$ value for each case.

###  IV.3 Bounded Cell

To study matrix-mediated mechanical interaction between cells in a nonlinear
medium, we employ the mean-field geometrical setup presented in Sec. II, which
we have studied for a linear matrix in Sec. III. We now numerically solve Eq.
(51) with the boundary condition $\tilde{u}(\tilde{R}_{1})=0$. In this bounded
geometry, the dimensionless displacement solution
$\tilde{u}\left(\tilde{r}\right)$ depends not only on the Poisson ratio $\nu$
and on the dimensionless nonlinearity $A$, but also on the dimensionless
distance $\tilde{R}_{1}$ between cells. We show in Fig. 3b numerical results
for varying $A$ and $\tilde{R}_{1}$.

The first thing to note in comparison to a linear medium is the following: In
the unbounded geometry, the dimensionless nonlinearity $A$ of the problem may
take arbitrarily large values. This means that for a given material
nonlinearity $b$, the displacement $u_{0}$ on the cell boundary can take any
value. As seen in Fig. 3a, increasing $A$ causes an increase in the radius
$\tilde{R}_{*}$ where $\tilde{u}\left(\tilde{r}\right)$ crosses over from the
compression-only solution $\tilde{u}=\tilde{r}$ to the shear-only solution
$\tilde{u}=\tilde{u}_{\rm eff}/\tilde{r}^{2}$. However, for the bounded
geometry, this crossover length $\tilde{R}_{*}$ clearly cannot be larger than
$\tilde{R}_{1}$. This implies that for given dimensionless cell-cell distance
$\tilde{R}_{1}$, there is a maximal value $A_{a}$ that the nonlinearity can
reach. Thus for given geometry specified by $R_{0}$ and $R_{1}$ and for given
material nonlinearity $b$, the displacement $u_{0}$ on the cell boundary can
increase only up to some maximal value
$u_{a}=R_{0}\sqrt{A_{a}(\tilde{R}_{1})/b}$.

In our nonlinear model (39), for large strain the resistance to shear becomes
increasingly larger than the resistance to compression. Hence, in mechanical
equilibrium, the shear component of the deformation tends to be minimal. For
our spherical geometry, requiring $\tilde{I}_{3}=0$ leads to the solution
$\tilde{u}\left(\tilde{r}\right)=\tilde{r}$, which is consistent with the
boundary condition $\tilde{u}(1)=1$, but cannot be consistent with the
boundary condition $\tilde{u}\left(\infty\right)=0$ or
$\tilde{u}(\tilde{R}_{1})=0$. An approximate solution to (51) that decays in
the far field is obtained by taking $\tilde{I}_{3}=1/A$ Shokef2012 , which by
substitution in (52) leads to

$\tilde{u}\left(\tilde{r}\right)=\tilde{r}-\sqrt{\frac{3}{4A}}\tilde{r}\log{\tilde{r}},$
(53)

which in turn also causes $d\tilde{I}_{3}/d\tilde{r}$ to vanish. Figure 3
shows how well (53) agrees with the near-field behavior of the displacement
field, both for the bounded and for the unbounded geometry.

Figure 4: (a) Maximal possible cell contraction increases with cell-cell
distance according to Eq. (56). Numerical results for
$b=0.1,0.5,0.8,1.1,2,5,10$ all collapse to a single curve, which perfectly
agrees with the theoretical prediction. (b) Position of maximal displacement
increases with dimensionless nonlinearity according to Eq. (59) both for the
bounded and for the unbounded geometries. (c) Effective far field displacement
increases with $A$ and decreases with $\tilde{R}_{1}$ following Eq. (61). In
all panels, theoretical prediction is marked by red dotted line.

We may approximate the maximal possible displacement $u_{a}$ on the cell
boundary by requiring that for $A_{a}$ the approximate near-field solution
(53) reaches zero at the distance $\tilde{R}_{1}$. The justification for this
is that for $u_{0}<u_{a}$, or correspondingly $A<A_{a}$, the displacement
field has room to cross over at larger distances to the linear elasticity
solution so that it will eventually satisfy the boundary condition
$\tilde{u}(\tilde{R}_{1})=0$. For $u_{0}>u_{a}$ on the other hand, the near-
field solution (53) is still finite at $\tilde{R}_{1}$ and the zero-
displacement boundary condition there may not be satisfied. Thus, setting
$\tilde{u}\left(\tilde{R}_{1}\right)=0$ in (53), we write

$\tilde{R}_{1}-\sqrt{\frac{3}{4A_{a}}}\tilde{R}_{1}\log{\tilde{R}_{1}}=0,$
(54)

from which we obtain the maximal dimensionless nonlinearity

$A_{a}=\frac{3}{4}\log^{2}\tilde{R_{1}}.$ (55)

This lead to the maximal possible cell contraction

$\tilde{u}_{a}=\frac{u_{a}}{R_{0}}=\sqrt{\frac{3}{4}}\gamma_{a}\log{\frac{R_{1}}{R_{0}}},$
(56)

which scales linearly with the maximal shear strain $\gamma_{a}$, and
logarithmically with the dimensionless cell-cell distance. This simple
theoretical prediction agrees remarkably well with our numerical results, as
shown in Fig. 4a. Note that at $A=A_{a}$, our approximate solution (53)
becomes exact, as it satisfies the boundary conditions both at $\tilde{r}=1$
and at $\tilde{r}=\tilde{R}_{1}$, see Fig. 3b. Before moving on to discussing
the interaction energy we also note that (56) may be inverted to obtain the
closest distance between cells that contract by a given amount $u_{0}$,

$\min(R_{1})=R_{0}\exp{\left(\sqrt{\frac{4}{3}}\frac{u_{0}/R_{0}}{\gamma_{a}}\right)},$
(57)

or alternatively

$\min\left(\tilde{R}_{1}\right)=\exp{\left(\sqrt{\frac{4A}{3}}\right)}.$ (58)

###  IV.4 Interaction Energy

To calculate the interaction energy for given cell size $R_{0}$, cell-cell
distance $R_{1}$, material nonlinearity $b$, and cell contraction $u_{0}$, we
need to first numerically solve the boundary value problem to get the
displacement field $u(r)$ for the bounded and for the unbounded geometries for
any intermediate cell contraction $0<w<u_{0}$. Then, we use $u$ and $du/dr$ at
$R_{0}$ to evaluate the radial stress (46) $\tau(w)=\sigma_{rr}(R_{0})$ on the
cell surface. Finally, we numerically integrate Eq. (1) to obtain the elastic
energy stored in the medium for the bounded and for the unbounded cases, and
subtract them to get the interaction energy.

Figure 5: Interaction energy vs. cell contraction for different values of
$R_{1}$ and $b$. Interaction energy is scaled by $\mu R_{0}$, and cell
displacement is scaled by $R_{0}$. Black dashed lines are linear elasticity
behavior for $b=0$.

####  IV.4.1 Displacement Regulation

The numerical results of $\Delta E$ vs $u_{0}$ for various values of $b$ and
$R_{1}$ are given in Fig. 5. For small $u_{0}$, the material responds
linearly, thus the interaction energy grows quadratically with $u_{0}$, as
predicted in a linear medium. For increasing $u_{0}$, the interaction energy
grows more rapidly and diverges at a finite contraction $u_{a}$, which depends
on $R_{1}$ and on $b$ according to (56). We first note that the dependence of
$\Delta E$ on the material nonlinearity $b$ may be scaled out if we normalize
$\Delta E$ by the interaction energy (15) that we would get for these values
of $u_{0}$, $R_{0}$ and $R_{1}$ in a linear medium. Figure 6 shows that after
this normalization, and by also scaling $u_{0}$ by its maximal possible value
$u_{a}$ in the nonlinear case, we obtain collapse of the results for different
$b$ values to master curves that depend only on $\tilde{R}_{1}=R_{1}/R_{0}$.

Figure 6: Scaled interaction energy behaves according to linear elasticity
(15) for small cell contraction (black dashed line). As $u_{0}$ approaches
$u_{a}$, the interaction energy grows according to Eq. (63) (red dashed
lines), until it diverges at $u_{0}=u_{a}$ (blue dashed line). Results for
different values of the material nonlinearity $b$ collapse to curves that
depend only on $\tilde{R}_{1}$.

To approximate the entire dependence on $u_{0}$ as well as the dependence on
$\tilde{R}_{1}$, we first note that the near-field approximate solution (53)
is maximal at

$\tilde{R}_{*}=\exp{\left(\sqrt{\frac{4A}{3}}-1\right)}.$ (59)

This agrees very well with the position of the maximum of
$\tilde{u}(\tilde{r})$ as obtained from our numerical solution, both for
bounded and for unbounded geometries, see Fig. 4b. Note that (53) has a local
maximum at $\tilde{r}>1$ only for $A>3/4$, thus all the analysis that follows
is relevant only for dimensionless nonlinearity that is larger than this
level. Now, in the far field, we expect the solution (12) in a linear medium
to hold, but with the displacement $u_{0}$ on the surface of the cell replaced
by some larger, effective displacement $u_{\rm eff}$. In dimensionless form
this read

$\tilde{u}\left(\tilde{r}\right)=\frac{\tilde{u}_{\rm
eff}}{\tilde{R}_{1}^{3}-1}\left(\frac{\tilde{R}_{1}^{3}}{\tilde{r}^{2}}-\tilde{r}\right).$
(60)

As has successfully been done for the unbounded geometry Shokef2012 , also for
our bounded geometry, the dimensionless effective displacement may be obtained
by matching the near field solution (53) and the far field solution (60) at
$\tilde{R}_{*}$ (59). The crossover between these two solutions is not sharp
but exists over a certain crossover region, see Fig. 7. Thus we allow for a
dimensionless multiplicative factor $\alpha$ between the two solutions at
$\tilde{R}_{*}$. This leads to

$\tilde{u}_{\rm eff}=\frac{u_{\rm
eff}}{u_{0}}=\alpha\sqrt{\frac{3}{4A}}\frac{\tilde{R}_{1}^{3}-1}{\frac{\tilde{R}_{1}^{3}}{\tilde{R}_{*}^{3}}-1}.$
(61)

For $\tilde{R}_{1}\gg\tilde{R}_{*}$ this reduces to $\tilde{u}_{\rm
eff}=\alpha\sqrt{3/(4A)}\tilde{R}_{*}^{3}$. Note that (61) cannot explain the
divergence of $\Delta E$ as $u_{0}$ approaches $u_{a}$ since it is not
singular at the maximal cell contraction $A_{a}$; it is just equal to

$\tilde{u}_{\rm
eff}=\alpha\frac{\tilde{R}_{1}^{3}-1}{\log\tilde{R}_{1}\left(e^{3}-1\right)}.$
(62)

Figure 4c shows the agreement of (61) with the numerical results, and
identifies the value $\alpha\approx 4.7$ for all the values of $b$ and
$\tilde{R}_{1}$ considered.

Figure 7: Crossover of $\tilde{u}(\tilde{r})$ from near-field behavior of Eq.
(53) (dashed black line) to far field $\tilde{u}=\tilde{u}_{\rm
eff}/\tilde{r}^{2}$ behavior. Matching the two solutions at $\tilde{R}_{*}$
(green dashed line) underestimates $\tilde{u}_{\rm eff}$. Matching at a
slightly larger distance (orange dashed line) is obtained by taking $\alpha>1$
in Eq. (61). Results are shown for A=35.7858.

The near field behavior of the displacement field is very similar in the
bounded and in the unbounded geometries, and most of the differences between
the two cases are in the far field (see Fig. 3). Thus we assume that most of
the contribution to the difference in the elastic energy stored in the medium
between the two geometries comes from the far field region. This enables us to
approximate the interaction energy as the difference in energies stored in a
linearly elastic medium displaced by $u_{\rm eff}$ rather than by $u_{0}$ at
the cell boundary. Thus we substitute $u_{\rm eff}$ of the bounded geometry in
Eq. (15) for $\Delta E$, and get

$\Delta E=\Delta E_{\rm
lin}\frac{3\alpha^{2}}{4A}\left(\frac{\tilde{R}_{1}^{3}-1}{\frac{\tilde{R}_{1}^{3}}{\tilde{R}_{*}^{3}}-1}\right)^{2},$
(63)

where $\Delta E_{\rm lin}$ is the interaction energy in a linear medium (15).
Figure 6 shows the impressive agreement of this theoretical expression with
the numerical results up until very close to the divergence of $\Delta E$ at
$u_{a}$.

####  IV.4.2 Stress Regulation

We now numerically evaluate the radial stress (46) on the surface of the cell
for the unbounded and for the bounded geometries, and subtract the
corresponding energies at given stress to obtain the interaction energy for
the situation in which the cell regulates the stress on its surface regardless
of its mechanical environment. For simplicity, for cells in a nonlinear matrix
we will assume that the cell does not resist the deformation and that all the
active force applied by the cell is directed toward generating stress in the
matrix. Namely $\sigma_{c}=0$ and $\sigma_{0}=-\sigma_{m}$. When relating this
to the linear elasticity solution (20), this is obtained by setting $K_{c}=0$.

The reasoning explaining why for displacement regulation the interaction
energy is positive (repulsion) and for stress regulation it is negative
(attraction) is valid also for cells in a strain stiffening matrix. Our
numerical calculations indeed result in positive interaction energy for
displacement regulation and negative interaction energy for stress regulation.
For displacement regulation we saw that the divergence of the shear stress at
a finite strain $\gamma_{a}$ (see Fig. 2a) causes the interaction energy to
diverge at a finite displacement $u_{a}$. For stress regulation on the other
hand, at least from the point of view of the matrix resistance, the active
stress that the cell generates can be arbitrarily large (see Fig. 2b).

Figure 8: Interaction energy with stress regulation in a nonlinear medium.
(a-c) Results for different values of $\tilde{R}_{1}$ and $b$. Dashed black
lines are linear elasticity behavior (20). Interaction energy is scaled by
$\mu/R_{0}^{3}$ and stress is scaled by $\mu$. (d) When normalizing the
interaction energy by the linear prediction and normalizing the stress by the
typical stress $\mu\gamma_{a}$ where strain stiffening becomes dominant,
results for all $b$ and $\tilde{R}_{1}$ collapse to a single curve.

When plotting in Fig. 8 the numerically obtained interaction energy vs. the
regulated stress we see the linear-elasticity $\Delta E\propto-\sigma_{0}^{2}$
behavior (20) at low stress. This is followed by a more complex crossover at
intermediate non-linearity, and eventually at extremely large stress the
interaction energy again scales as $\sigma_{0}^{2}$, enabling arbitrarily
large values of the stress, as expected. Interestingly, we obtain the
quadratic scaling $\Delta E\propto\sigma_{0}^{2}$ not only in the linear
regime but also in the strongly nonlinear regime. This may be explained from
the fact that the nonlinear dependence of $u_{\rm eff}$ on $u_{0}$ and of
$\sigma_{0}$ on $u_{\rm eff}$ are reciprocal such that $\sigma_{0}\propto
u_{\rm eff}$ Shokef2012 . Combining this with the result presented above that
$\Delta E\propto u_{\rm eff}^{2}$ leads to $\Delta E\propto\sigma_{0}^{2}$.
Impressively, when we scale the nonlinear interaction energy by the linear
result, and scale the stress by the typical stress $\mu\gamma_{a}$ where the
strain stiffening becomes significant, all the numerical results collapse to a
single curve, as shown in Fig. 8d for multiple values of $\tilde{R}_{1}$ and
$b$.

##  V Discussion

In this paper, we used a mean-field approach to study matrix mediated
interactions between contracting biological cells in linear and in nonlinear
elastic surroundings. We showed how the regulatory behavior of the cell’s
mechanical activity always gives rise to repulsion when displacement is
regulated and to attraction when stress is regulated. For displacement
regulation, the interaction does not depend on the rigidity of the cell. For
stress regulation, it depends only on the bulk modulus of the cell in
spherically symmetric situations, while in cylindrical setups it depends also
on its shear modulus. For a nonlinear, shear stiffening matrix, the
interaction energy diverges at a finite displacement, while for stress
regulation the stress may be arbitrarily large.

For simplicity, we treated the cell as a passive, linearly elastic solid. This
is clearly a very crude description, and the way that cells set their size,
displacement and forces include many further processes, see e.g. Adar2020 ;
Doss2020 . However, we emphasize that for displacement regulation the
mechanical response of the cell is irrelevant for the interaction energy.
Namely, no matter how complicated is the response of the cell, if we consider
the situation in which cells regulate the displacement they generate, then the
inside of the cell behaves exactly the same regardless of the distance to
neighboring cells, and our results both for a linearly elastic matrix and for
a nonlinear matrix are valid also beyond the simplified assumption that the
cell is a passive, linearly elastic solid. In the case of stress regulation,
our analytical results for a linearly elastic matrix rely on the linearity
assumption for the cell, whereas for a nonlinear matrix we assumed that the
cell is much softer than the matrix, and therefore its mechanical response
does not affect the interaction with other cells.

We focused our nonlinear material model for the extracellular matrix on shear
stiffening and included nonlinearity only in the resistance to shear. It would
be interesting to study the interactions and different regulatory behaviors
that we consider also in models that take into account the anisotropy of the
matrix Xu2015 ; Wang2020 ; Goren2020 , that include also nonlinearity in the
compressive response, and in models that take into account the discrete
fibrous nature of the biopolymer gel comprising the extracellular environment
Sopher2018 ; Goren2020 ; Natan2021 ; Alisafaei2021 .

###### Acknowledgements.

We thank Dan Ben-Yaakov, Roman Golkov, Shahar Goren, Yoni Koren, Ayelet Lesman
and Sam Safran for fruitful discussions. This research was supported in part
by a grant from the United States-Israel Binational Science Foundation, by the
Israel Science Foundation Grant No. 968/16 and by the National Science
Foundation Grant No. NSF PHY-1748958.

## References

  * (1) U. S. Schwarz and S. A. Safran, Rev. Mod. Phys. 85, 1327 (2013).
  * (2) A. Engler, S. Sen, H. Sweeney, and D. Discher, Cell 126, 677 (2006).
  * (3) M. Poujade, E. Grasland-Mongrain, A. Hertzog, J. Jouanneau, P. Chavrier, B. Ladoux, A. Buguin, and P. Silberzan, Proc. Natl. Acad. Sci. USA 104, 15988 (2007).
  * (4) D. Montell, Nat. Rev. Mol. Cell Biol. 4, 13 (2003).
  * (5) A. Lesman, J. Notbohm, D. A. Tirrell, and G. Ravichandran, J. Cell Biol. 205, 155 (2014).
  * (6) A. Abuhattum, A. Gefen, and D. Weihs, Integr. Biol. 7, 1212 (2015).
  * (7) P. Friedl and D. Gilmour, Nat. Rev. Mol. Cell Biol. 10, 445 (2009).
  * (8) N. Gal and D. Weihs, Cell Biochem. Biophys. 63, 199 (2012).
  * (9) D. Gomez, S. Natan, Y. Shokef, and A. Lesman, Adv. Biosys. 3, 1900192 (2019).
  * (10) D. Gomez, E. Teomy, A. Lesman, and Y. Shokef, New J. Phys. 22, 103008 (2020).
  * (11) W-H. Jung, N. Yam, C-C. Chen, K. Elawad, B. Hu, and Y. Chen, Biomaterials 234, 119756 (2020).
  * (12) C. A. Reinhart-King, M. Dembo, and D. A. Hammer, Biophys. J. 95, 6044 (2008).
  * (13) J. P. Winer .S. Oake, and P. A. Janmey, PLoS ONE 4, e6382 (2009).
  * (14) A. Buxboim, K. Rajagopal, A. E. X. Brown, and D. E. Discher, J. Phys.: Condens. Matter 22, 194116 (2010).
  * (15) H. Mohammadi, P. A. Janmey, C. A. McCulloch, Biomaterials 35, 1138 (2014).
  * (16) J. Notbohm, A. Lesman, P. Rosakis, D. A. Tirrell, and G. Ravichandran, J. R. Soc. Interface 12, 20150320 (2015).
  * (17) J. Notbohm, A. Lesman, D. A. Tirrell, and G. Ravichandran, Integr. Biol. 7, 1186 (2015).
  * (18) U. S. Schwarz and S. A. Safran, Phys. Rev. Lett. 88, 048102 (2002).
  * (19) K. Lau and W. Kohn, Surface Science 65, 607 (1977).
  * (20) I. B. Bischofs and U. S. Schwarz, Proc. Natl. Acad. Sci. USA 100, 9274 (2003).
  * (21) I. B. Bischofs, S. A. Safran, and U. S. Schwarz, Phys. Rev. E 69, 021911 (2004).
  * (22) I. B. Bischofs and U. S. Schwarz, Phys. Rev. Lett. 95, 068102 (2005).
  * (23) I. B. Bischofs and U. S. Schwarz, Acta Biomater. 2, 0510391 (2006).
  * (24) R. Golkov and Y. Shokef, Phys. Rev. E 99, 032418 (2019).
  * (25) Y. Shokef and S. A. Safran, Phys. Rev. Lett. 108, 178103 (2012); Phys. Rev. Lett. 109, 169901 (2012).
  * (26) X. Xu and S. A. Safran, Phys. Rev. E 92, 032728 (2015).
  * (27) D. Ben-Yaakov, R. Golkov, Y. Shokef, and S. A. Safran, Soft Matter 11, 1412 (2015).
  * (28) H. Wang and X. Xu, Soft Matter 16, 10781 (2020).
  * (29) J. N. Israelachvili, Intermolecular and Surface Forces (Academic Press, London, 1992)
  * (30) R. Golkov and Y. Shokef, New J. Phys. 19, 063011 (2017).
  * (31) R. S. Sopher, H. Tokash, S. Natan, M. Sharabi, O. Shelah, O. Tchaicheeyan, and A. Lesman, Biophys. J. 115, 1357 (2018).
  * (32) R. Brown, R. Prajapati, D. McGrouther, I. Yannas, and M. Eastwood, J. Cell. Physiol. 175, 323 (1998).
  * (33) T. M. Freyman, I. V. Yannas, R. Yokoo, and L. J. Gibson, Exp. Cell. Res. 272, 153 (2002).
  * (34) A. Saez, A. Buguin, P. Silberzan, and B. Ladoux, Biophys. J. 89, L52 (2005).
  * (35) M. Ghibaudo, A. Saez, L. Trichet, A. Xayaphoummine, J. Browaeys, P. Silberzan, A. Buguin, and B. Ladoux, Soft Matter 4, 1836 (2008).
  * (36) R. De, A. Zemel, and S. A. Safran, Biophys. J. 94, L29 (2008).
  * (37) S. He, Y. Su, B. Ji, and H. Gao, J. Mech. Phys. Solids 70, 116 (2014).
  * (38) S. Natan, Y. Koren, O. Shelah, S. Goren, and A. Lesman, Mol. Biol. Cell 31, 1474 (2020).
  * (39) A. Bower, Applied Mechanics of Solids (CCRC Press, Boca Raton, FL, USA, 2009).
  * (40) M. L. Gardel, J. H. Shin, F. C. MacKintosh, L. Mahadevan, P. Matsudaira, and D. A. Weitz, Science 304, 1301 (2004).
  * (41) C. Storm, J. J. Pastore, F. C. MacKintosh, T. C. Lubensky, and P. A. Janmey, Nature 435, 191 (2005).
  * (42) D. Vader, A. Kabla, D. Weitz, L. Mahadevan, PLoS ONE 4, e5902 (2009).
  * (43) C. P. Broedersz and F. C. MacKintosh, Rev. Mod. Phys. 86, 995 (2014).
  * (44) R. M. Adar and S. A. Safran, Proc. Natl. Acad. Sci. USA 117, 5604 (2020).
  * (45) B. L. Doss, M. Pan, M. Gupta, G. Grenci, R-M. Mège, C. T. Lim, M. P. Sheetz, R. Voituriez, and B. Ladoux, Proc. Natl. Acad. Sci. USA 117, 12817 (2020).
  * (46) S. Goren, Y. Koren, X. Xu, and A. Lesman, Biophys. J. 118, 1152 (2020).
  * (47) F. Alisafaei, X. Chen, T. Leahy, P. A. Janmey, and V. B. Shenoy, Soft Matter 17, 241 (2021).


employing the IRAC $8~{}\micron$, MIPS $24~{}\micron$, PACS $70~{}\micron$,
and SPIRE $500~{}\micron$ bands. This is inspired by the work carried out by
Calzetti et al. (2018); Gregg et al. (2022), where this colour–colour plot is
generated at a kpc resolution for nearby galaxies. For the rest of this
section, we use the notation $F_{n}$ to represent $\nu F_{\nu}(n)$, the flux
in a band with pivot wavelength $n$.

Fig. 21 shows how the two libraries cover this plane. We also show the
relation found by Calzetti et al. (2018) along with the colors for individual
pixels for the “Regime-1” galaxies discussed in Gregg et al. (2022) which
exhibit uniform and high star-formation surface densities. Note that the curve
from Calzetti et al. (2018) uses $F_{1100}$ instead of $F_{500}$, here we have
made the conversion based on the factor given in Gregg et al. (2022).

As we compare the trends on this color–color plot, it is worth noting that
$F_{500}$ and $F_{8}$ are impacted by the dust emission originating outside of
the star-forming regions. The observational sample, despite being the one with
high star formation density is affected by the presence of cold dust
($F_{500}$) and to some extent, dust heating by old stars ($F_{8}$). This
model data will move downwards and slightly rightwards once the dust emission
originating outside of our model star-formation regions is taken into account.
This trend is also expected based on the nature of the Calzetti et al. (2018)
curve, where the addition of diffuse, cooler dust moves the data points
rightwards and downwards.

In the models considered in Fig. 21, an increase in $n_{\mathrm{cl}}$ and/or
$\epsilon_{\mathrm{SF}}$ leads to an increase in the luminosity-weighted
$\mathcal{F}_{\mathrm{in}}$ resulting in an increase in $F_{70}/F_{500}$, a
trend observed with increasing $\log C$ (and to some extent with increasing
$f_{\mathrm{PDR}}$) in the case of HiiM3. The lower $\epsilon_{\mathrm{SF}}$
systems tend to have a higher number of shell re-collapse events leading to an
increasingly narrow range of values in $F_{70}/F_{500}$ as the effective star
formation efficiency ($M_{\star}/M_{\mathrm{cl}}$ at the end of 30 Myr) could
reach values closer to higher $\epsilon_{\mathrm{SF}}$ models. The HiiM3
parameter, $f_{\mathrm{PDR}}$, is completely free from the evolution of the
bubbles and serves as a geometrical factor on the non-ionizing UV absorbed by
the molecular cloud assuming a fixed column depth. At a fixed $\log C$, this
leads to models moving rightward and downward on the MIR–FIR plane due to the
increasing contribution of PAHs and cold dust. In contrast, TODDLERS calculate
IR emission in a single model, and column depths are affected by the state of
the component shells of the star-forming complexes. The youngest and the most
massive components tend to have the highest impact on the IR SED. These
components are ones that are still embedded in their birth cloud and have high
dust column depths. If we consider $F_{8}/F_{24}$ a proxy for
$q_{\mathrm{PAH}}$, at lower densities, the unswept clouds are more diffused
and contain neutral Hydrogen. This leads to a higher luminosity-weighted PAH
abundance, moving the models rightward with decreasing $n_{\mathrm{cl}}$ at a
fixed $\epsilon_{\mathrm{SF}}$. Increasing $\epsilon_{\mathrm{SF}}$ at a fixed
$n_{\mathrm{cl}}$ tends to decrease the PAH fraction, leading to leftward
movement on the MIR–FIR plane.

On this MIR–FIR color plane the two libraries cover markedly different
regions, with the TODDLERS library showing a closer match to the observational
colors. We confirm that this trend is also seen in the colors generated with
simulated galaxies accounting for dust outside the star-forming regions
(subject of Paper 2 in this series). We attribute the large differences in the
$F_{8}/F_{24}$ to the differences in the dust models used in the two
libraries, in particular, the dust model employed in HiiM3 appears to be more
emissive around the $24\,\micron$ band.

##  7 Summary and outlook

In this work, we have presented a new emission library TODDLERS with the
primary aim of producing time-dependent emission diagnostics from gas and dust
around young stars for simulated galaxies. To achieve this, we have run a
large suite of semi-analytic calculations that allow us to infer the gas–star
geometry as the gas evolves under stellar feedback. The calculations assume a
cloud with a constant density profile and a finite mass evolving under the
influence of the feedback of a central cluster. This idealized approach allows
us to sweep a large parameter space while accounting for complex feedback
physics in a simplified fashion. The dominant stellar feedback channel evolves
as a function of metallicity. At high metallicities, the gas is pushed
predominantly by stellar winds and the subsequent SNe, whereas, the dominant
feedback channel at metallicities below $Z\lesssim 0.004$ is the Ly$\alpha$
radiation pressure due to multiple resonant scatterings. This is especially
the case for the higher end of cloud densities and masses where there is ample
neutral gas for prolonged periods. At higher metallicities, the role of
Ly$\alpha$ radiation pressure is subdominant but non-negligible. The clues
from this idealized Ly$\alpha$ radiation pressure setup point to the
importance of this feedback channel in young star-forming clouds warranting
more detailed studies, even at higher metallicities.

The semi-analytic calculations are passed on to Cloudy to carry out
calculations involving detailed chemistry enabling us to produce time-
dependent UV–mm observables. The tabulated observables include emission lines
originating from ionized, photo-dissociation, and molecular regions along with
the nebular, stellar, and dust continuum emission. We have used the BPT
diagram and the IRAS color-color diagram to map the parameter space of the
models onto the observable space demonstrating that they populate the expected
regions when compared to observational data from nearby galaxies and the Milky
Way. We have integrated these observables into SKIRT to use these data for
post-processing simulations where star-forming regions are not resolved by
assuming a power-law cloud population. A comparison focused on the IR colors
produced using the aforementioned TODDLERS implementation to that of the
currently only available star-forming regions’ library in SKIRT was carried
out. When confronted with observations, TODDLERS’ IR colors are in better
agreement with observations in comparison to HiiM3, which until now was the
only option in SKIRT to incorporate star forming regions’ emission. In a
companion paper, we use TODDLERS to produce UV-submm diagnostics using
simulated galaxies. By doing this, we can make a direct and detailed
comparison between broadband and line-emission data from simulated galaxies
with those from recent observational studies, such as Gregg et al. (2022) and
Groves et al. (2023).

This work serves as a proof of concept where we followed the evolution of
homogeneous clouds using a single set of stellar templates and the observable
generation employs a single chemical abundance set and dust model. Thus,
several areas remain open for exploration within the framework employed in
this work. As far as the evolutionary model is concerned, the cloud density
profile dictates its binding energy, modifying which could lead to significant
changes in the system’s evolution (Rahner et al., 2019). For example, consider
two clouds with identical mean density and mass. One exhibits a Bonnor-Ebert
profile with a high-density core (Ebert, 1955; Bonnor, 1956), while the other
has a constant density. Comparatively, the cloud with the Bonnor-Ebert profile
possesses a greater gravitational binding energy than its constant-density
counterpart, making it more resistant to destruction.

We have assumed that the unswept cloud is not dynamically affected by stellar
feedback, which, for example, lowers the effective outward momentum deposition
at low metallicities when the ionization front lies outside the shell. We have
also neglected the density and velocity gradients that would appear as a
result of the Ly$\alpha$ feedback which could lower the gas columns seen by
Ly$\alpha$ photons in non-trivial ways. While the problem of how feedback
affects the gas in star-forming regions is intrinsically 3D, relaxing these
assumptions in spherical symmetry would be an interesting iteration of the
current work (see, for example, Kourniotis et al., 2023).

We have shown that one of the key parameters that determine the emission
properties of a nebula, $U$, is critically dependent on the properties
resulting directly from the stellar templates ($Q_{\mathrm{H}}$ and
$F_{\mathrm{ram}}$). At the same time, the dust temperature is a function of
$\mathcal{F}_{\mathrm{in}}$, which depends on the stellar feedback. Such
dependencies demand a thorough investigation of the results presented here
while changing the stellar library. Modifying the stellar library could mean,
among other things, varying the initial elemental abundances, the physical
processes employed in the stellar models, and its IMF. Grasha et al. (2021)
have shown the need to use abundance sets beyond the metallicity-scaled solar
abundances for massive stars, especially due to their effects on the line
emission at the lower metallicity end. Similarly, the presence of stellar
rotation, binarity, and changes in the IMF play important roles in controlling
the ionizing photon production rate and mass loss from the stellar clusters
(Leitherer et al., 2014; Stanway & Eldridge, 2019). Such model choices will
have an impact on the various feedback channels incorporated here, leading to
interesting effects on the emission lines and dust emission. Thus considering
variations of the stellar templates based on the abundance sets, physical
processes, and IMF represent interesting avenues to explore.

Another intriguing prospect is to consider variations in the dust-to-metal
ratio as a function of metallicity and environment, which are expected based
on numerical modeling and observations (Asano et al., 2013; Rémy-Ruyer et al.,
2013; Schneider et al., 2016; De Vis et al., 2019). The computational
efficiency offered by the 1D calculations makes it feasible to generate models
with the aforementioned variations, we plan to pursue this in the future.

##  8 Data Availability

The TODDLERS library for post-processing galaxy formation simulations is
available for download along with the SKIRT code. All other data used in this
work, including the Cloudy output for individual models, are publicly
available at www.toddlers.ugent.be.

##  9 Acknowledgements

We thank Ilse De Looze, Jérémy Chastenet, and Brian Van Den Noortgaete for
useful discussions. We thank Benjamin Gregg for providing the MIR-FIR colors’
data from the KINGFISH sample. We thank Qing-Zeng Yan for providing the IRAS
colors’ data. AUK, AN, MB acknowledge the financial support of the Flemish
Fund for Scientific Research (FWO-Vlaanderen). The simulations carried out for
this work used the Tier-2 facilities of the Flemish Supercomputer Center
(https://www.vscentrum.be/) located at Ghent University. We are grateful to
the members, particularly Gary Ferland and Peter van Hoof, for addressing
numerous questions on the Cloudy online forum. We also thank the anonymous
referee for their valuable comments and suggestions.

## References

  * Abel et al. (2008) Abel N. P., van Hoof P. A. M., Shaw G., Ferland G. J., Elwert T., 2008, ApJ, 686, 1125
  * Asano et al. (2013) Asano R. S., Takeuchi T. T., Hirashita H., Inoue A. K., 2013, Earth, Planets and Space, 65, 213
  * Baes et al. (2011) Baes M., Verstappen J., De Looze I., Fritz J., Saftly W., Vidal Pérez E., Stalevski M., Valcke S., 2011, ApJS, 196, 22
  * Baes et al. (2019) Baes M., Trčka A., Camps P., Nersesian A., Trayford J., Theuns T., Dobbels W., 2019, MNRAS, 484, 4069
  * Baldwin et al. (1981) Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5
  * Barrientos Acevedo et al. (2023) Barrientos Acevedo D., et al., 2023, MNRAS, 524, 907
  * Belfiore et al. (2016) Belfiore F., et al., 2016, MNRAS, 461, 3111
  * Bonnor (1956) Bonnor W. B., 1956, MNRAS, 116, 351
  * Boquien et al. (2019) Boquien M., Burgarella D., Roehlly Y., Buat V., Ciesla L., Corre D., Inoue A. K., Salas H., 2019, A&A, 622, A103
  * Burkhart et al. (2015) Burkhart B., Lee M.-Y., Murray C. E., Stanimirović S., 2015, ApJ, 811, L28
  * Byler et al. (2017) Byler N., Dalcanton J. J., Conroy C., Johnson B. D., 2017, ApJ, 840, 44
  * Calzetti et al. (2018) Calzetti D., et al., 2018, ApJ, 852, 106
  * Camps & Baes (2015) Camps P., Baes M., 2015, Astronomy and Computing, 9, 20
  * Camps & Baes (2020) Camps P., Baes M., 2020, Astronomy and Computing, 31, 100381
  * Camps et al. (2016) Camps P., Trayford J. W., Baes M., Theuns T., Schaller M., Schaye J., 2016, MNRAS, 462, 1057
  * Camps et al. (2021) Camps P., Behrens C., Baes M., Kapoor A. U., Grand R., 2021, ApJ, 916, 39
  * Camps et al. (2022) Camps P., Kapoor A. U., Trcka A., Font A. S., McCarthy I. G., Trayford J., Baes M., 2022, MNRAS, 512, 2728
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chevalier & Klein (1978) Chevalier R. A., Klein R. I., 1978, ApJ, 219, 994
  * Chiang et al. (2018) Chiang I.-D., Sandstrom K. M., Chastenet J., Johnson L. C., Leroy A. K., Utomo D., 2018, ApJ, 865, 117
  * Churchwell (2002) Churchwell E., 2002, ARA&A, 40, 27
  * Churchwell et al. (2009) Churchwell E., et al., 2009, PASP, 121, 213
  * Colombo et al. (2014) Colombo D., et al., 2014, ApJ, 784, 3
  * Conroy (2013) Conroy C., 2013, ARA&A, 51, 393
  * De Vis et al. (2019) De Vis P., et al., 2019, A&A, 623, A5
  * Diemer et al. (2019) Diemer B., et al., 2019, MNRAS, 487, 1529
  * Dijkstra & Loeb (2008) Dijkstra M., Loeb A., 2008, MNRAS, 391, 457
  * Dijkstra & Loeb (2009) Dijkstra M., Loeb A., 2009, MNRAS, 396, 377
  * Dopita et al. (2002) Dopita M. A., Groves B. A., Sutherland R. S., Binette L., Cecil G., 2002, ApJ, 572, 753
  * Dopita et al. (2005) Dopita M. A., et al., 2005, ApJ, 619, 755
  * Dopita et al. (2006) Dopita M. A., et al., 2006, ApJ, 647, 244
  * Dopita et al. (2013) Dopita M. A., Sutherland R. S., Nicholls D. C., Kewley L. J., Vogt F. P. A., 2013, ApJS, 208, 10
  * Draine (1978) Draine B. T., 1978, ApJS, 36, 595
  * Draine (2011) Draine B. T., 2011, ApJ, 732, 100
  * Duffell (2016) Duffell P. C., 2016, ApJ, 821, 76
  * Ebert (1955) Ebert R., 1955, Z. Astrophys., 37, 217
  * El-Badry et al. (2019) El-Badry K., Ostriker E. C., Kim C.-G., Quataert E., Weisz D. R., 2019, MNRAS, 490, 1961
  * Elmegreen (2011) Elmegreen B. G., 2011, in Charbonnel C., Montmerle T., eds, EAS Publications Series Vol. 51, EAS Publications Series. pp 45–58 (arXiv:1101.3112), doi:10.1051/eas/1151004
  * Feldmann et al. (2022) Feldmann R., et al., 2022, arXiv e-prints, p. arXiv:2205.15325
  * Ferland et al. (2017) Ferland G. J., et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385
  * Fielding et al. (2020) Fielding D. B., Ostriker E. C., Bryan G. L., Jermyn A. S., 2020, ApJ, 894, L24
  * Franco et al. (1994) Franco J., Shore S. N., Tenorio-Tagle G., 1994, ApJ, 436, 795
  * Gebek et al. (2023) Gebek A., et al., 2023, MNRAS,
  * Geen & de Koter (2022) Geen S., de Koter A., 2022, MNRAS, 509, 4498
  * Grand et al. (2017) Grand R. J. J., et al., 2017, MNRAS, 467, 179
  * Grasha et al. (2021) Grasha K., Roy A., Sutherland R. S., Kewley L. J., 2021, ApJ, 908, 241
  * Grasha et al. (2022) Grasha K., et al., 2022, ApJ, 929, 118
  * Gregg et al. (2022) Gregg B., Calzetti D., Heyer M., 2022, ApJ, 928, 120
  * Grevesse et al. (2010) Grevesse N., Asplund M., Sauval A. J., Scott P., 2010, Ap&SS, 328, 179
  * Gronke et al. (2017) Gronke M., Dijkstra M., McCourt M., Oh S. P., 2017, A&A, 607, A71
  * Groves et al. (2008) Groves B., Dopita M. A., Sutherland R. S., Kewley L. J., Fischera J., Leitherer C., Brandl B., van Breugel W., 2008, ApJS, 176, 438
  * Groves et al. (2023) Groves B., et al., 2023, MNRAS, 520, 4902
  * Guidi et al. (2015) Guidi G., Scannapieco C., Walcher C. J., 2015, MNRAS, 454, 2381
  * Hanaoka et al. (2019) Hanaoka M., et al., 2019, PASJ, 71, 6
  * Harper-Clark & Murray (2009) Harper-Clark E., Murray N., 2009, ApJ, 693, 1696
  * Heyer & Dame (2015) Heyer M., Dame T. M., 2015, ARA&A, 53, 583
  * Heyer et al. (2009) Heyer M., Krawczyk C., Duval J., Jackson J. M., 2009, ApJ, 699, 1092
  * Hillier & Miller (1998) Hillier D. J., Miller D. L., 1998, ApJ, 496, 407
  * Hirschmann et al. (2022) Hirschmann M., et al., 2022, arXiv e-prints, p. arXiv:2212.02522
  * Hopkins et al. (2023) Hopkins P. F., et al., 2023, MNRAS, 519, 3154
  * Hunt & Hirashita (2009) Hunt L. K., Hirashita H., 2009, A&A, 507, 1327
  * Inoue (2002) Inoue A. K., 2002, ApJ, 570, 688
  * Jang et al. (2022) Jang J. K., et al., 2022, arXiv e-prints, p. arXiv:2211.00931
  * Jonsson et al. (2010) Jonsson P., Groves B. A., Cox T. J., 2010, MNRAS, 403, 17
  * Kannan et al. (2020) Kannan R., Marinacci F., Vogelsberger M., Sales L. V., Torrey P., Springel V., Hernquist L., 2020, MNRAS, 499, 5732
  * Kapoor et al. (2021) Kapoor A. U., et al., 2021, MNRAS, 506, 5703
  * Katz et al. (2022) Katz H., et al., 2022, arXiv e-prints, p. arXiv:2211.04626
  * Kauffmann et al. (2003) Kauffmann G., et al., 2003, MNRAS, 346, 1055
  * Kewley & Dopita (2002) Kewley L. J., Dopita M. A., 2002, ApJS, 142, 35
  * Kewley & Ellison (2008) Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183
  * Kewley et al. (2001) Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J., 2001, ApJ, 556, 121
  * Kewley et al. (2019) Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511
  * Kim et al. (2016) Kim J.-G., Kim W.-T., Ostriker E. C., 2016, ApJ, 819, 137
  * Kimm et al. (2018) Kimm T., Haehnelt M., Blaizot J., Katz H., Michel-Dansac L., Garel T., Rosdahl J., Teyssier R., 2018, MNRAS, 475, 4617
  * Kobayashi et al. (2022) Kobayashi M. I. N., Inoue T., Tomida K., Iwasaki K., Nakatsugawa H., 2022, ApJ, 930, 76
  * Kourniotis et al. (2023) Kourniotis M., Wünsch R., Martínez-González S., Palouš J., Tenorio-Tagle G., Ehlerová S., 2023, MNRAS, 521, 5686
  * Kreckel et al. (2019) Kreckel K., et al., 2019, ApJ, 887, 80
  * Krumholz (2013) Krumholz M. R., 2013, MNRAS, 436, 2747
  * Lancaster et al. (2021a) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021a, ApJ, 914, 89
  * Lancaster et al. (2021b) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021b, ApJ, 914, 90
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJS, 123, 3
  * Leitherer et al. (2014) Leitherer C., Ekström S., Meynet G., Schaerer D., Agienko K. B., Levesque E. M., 2014, ApJS, 212, 14
  * Leja et al. (2017) Leja J., Johnson B. D., Conroy C., van Dokkum P. G., Byler N., 2017, ApJ, 837, 170
  * Leja et al. (2019) Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019, ApJ, 876, 3
  * Levesque et al. (2010) Levesque E. M., Kewley L. J., Larson K. L., 2010, AJ, 139, 712
  * Li & Draine (2001) Li A., Draine B. T., 2001, ApJ, 554, 778
  * Lopez et al. (2011) Lopez L. A., Krumholz M. R., Bolatto A. D., Prochaska J. X., Ramirez-Ruiz E., 2011, ApJ, 731, 91
  * Maeder & Conti (1994) Maeder A., Conti P. S., 1994, ARA&A, 32, 227
  * Martínez-González et al. (2014) Martínez-González S., Silich S., Tenorio-Tagle G., 2014, ApJ, 785, 164
  * Matsumoto et al. (2023) Matsumoto K., Camps P., Baes M., De Ceuster F., Wada K., Nakagawa T., Nagamine K., 2023, Self-consistent dust and non-LTE line radiative transfer with SKIRT
  * Meynet et al. (1994) Meynet G., Maeder A., Schaller G., Schaerer D., Charbonnel C., 1994, A&AS, 103, 97
  * Miura et al. (2012) Miura R. E., et al., 2012, ApJ, 761, 37
  * Miville-Deschênes et al. (2017) Miville-Deschênes M.-A., Murray N., Lee E. J., 2017, ApJ, 834, 57
  * Moy et al. (2001) Moy E., Rocca-Volmerange B., Fioc M., 2001, A&A, 365, 347
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei. University Science Books
  * Ostriker & Cowie (1981) Ostriker J. P., Cowie L. L., 1981, ApJ, 243, L127
  * Pauldrach et al. (2001) Pauldrach A. W. A., Hoffmann T. L., Lennon M., 2001, A&A, 375, 161
  * Pellegrini et al. (2020) Pellegrini E. W., Rahner D., Reissl S., Glover S. C. O., Klessen R. S., Rousseau-Nepton L., Herrera-Camus R., 2020, MNRAS, 496, 339
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Pilyugin & Thuan (2005) Pilyugin L. S., Thuan T. X., 2005, ApJ, 631, 231
  * Popping et al. (2017) Popping G., Puglisi A., Norman C. A., 2017, MNRAS, 472, 2315
  * Popping et al. (2021) Popping G., et al., 2021, The dust-continuum size of TNG50 galaxies at $z=1-5$: a comparison with the distribution of stellar light, stars, dust and H2 (arXiv:2101.12218)
  * Priestley et al. (2022) Priestley F. D., De Looze I., Barlow M. J., 2022, MNRAS, 509, L6
  * Rahner et al. (2017) Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2017, MNRAS, 470, 4453
  * Rahner et al. (2018) Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2018, MNRAS, 473, L11
  * Rahner et al. (2019) Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2019, MNRAS, 483, 2547
  * Rémy-Ruyer et al. (2013) Rémy-Ruyer A., et al., 2013, A&A, 557, A95
  * Rémy-Ruyer et al. (2014) Rémy-Ruyer A., et al., 2014, A&A, 563, A31
  * Robitaille et al. (2012) Robitaille T. P., Churchwell E., Benjamin R. A., Whitney B. A., Wood K., Babler B. L., Meade M. R., 2012, A&A, 545, A39
  * Rodriguez-Gomez et al. (2019) Rodriguez-Gomez V., et al., 2019, MNRAS, 483, 4140
  * Rousseau-Nepton et al. (2018) Rousseau-Nepton L., Robert C., Martin R. P., Drissen L., Martin T., 2018, MNRAS, 477, 4152
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schneider et al. (2016) Schneider R., Hunt L., Valiante R., 2016, MNRAS, 457, 1842
  * Semenov et al. (2003) Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611
  * Skinner & Ostriker (2015) Skinner M. A., Ostriker E. C., 2015, ApJ, 809, 187
  * Smith & Hayward (2018) Smith D. J. B., Hayward C. C., 2018, MNRAS, 476, 1705
  * Smith et al. (2017) Smith A., Bromm V., Loeb A., 2017, MNRAS, 464, 2963
  * Smith et al. (2019) Smith A., Ma X., Bromm V., Finkelstein S. L., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2019, MNRAS, 484, 39
  * Smith et al. (2020) Smith A., Kannan R., Tsang B. T. H., Vogelsberger M., Pakmor R., 2020, ApJ, 905, 27
  * Smith et al. (2022) Smith A., et al., 2022, MNRAS, 517, 1
  * Smith et al. (2023) Smith M. C., et al., 2023, arXiv e-prints, p. arXiv:2301.07116
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Stanway & Eldridge (2019) Stanway E. R., Eldridge J. J., 2019, A&A, 621, A105
  * Tacchella et al. (2022) Tacchella S., et al., 2022, MNRAS, 513, 2904
  * Tan et al. (2021) Tan B., Oh S. P., Gronke M., 2021, MNRAS, 502, 3179
  * Tomaselli & Ferrara (2021) Tomaselli G. M., Ferrara A., 2021, MNRAS, 504, 89
  * Torrey et al. (2015) Torrey P., et al., 2015, MNRAS, 447, 2753
  * Tress et al. (2020) Tress R. G., Smith R. J., Sormani M. C., Glover S. C. O., Klessen R. S., Mac Low M.-M., Clark P. C., 2020, MNRAS, 492, 2973
  * Trčka et al. (2020) Trčka A., et al., 2020, MNRAS, 494, 2823
  * Trčka et al. (2022) Trčka A., et al., 2022, MNRAS, 516, 3728
  * Utomo et al. (2018) Utomo D., et al., 2018, ApJ, 861, L18
  * Vale Asari et al. (2020) Vale Asari N., et al., 2020, MNRAS, 498, 4205
  * Vander Meulen et al. (2023) Vander Meulen B., Camps P., Stalevski M., Baes M., 2023, arXiv e-prints, p. arXiv:2304.10563
  * Verdolini et al. (2013) Verdolini S., Yeh S. C. C., Krumholz M. R., Matzner C. D., Tielens A. G. G. M., 2013, ApJ, 769, 12
  * Vogelsberger et al. (2020a) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020a, Nature Reviews Physics, 2, 42
  * Vogelsberger et al. (2020b) Vogelsberger M., et al., 2020b, MNRAS, 492, 5167
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Weingartner & Draine (2001) Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296
  * Yan et al. (2018) Yan Q.-Z., et al., 2018, MNRAS, 476, 3981
  * Yang et al. (2023) Yang S., Lidz A., Smith A., Benson A., Li H., 2023, arXiv e-prints, p. arXiv:2304.09261
  * Zavagno et al. (2010) Zavagno A., et al., 2010, A&A, 518, L101
  * van Hoof (2022) van Hoof P. A. M., 2022, Private Communication
  * van Hoof et al. (2004) van Hoof P. A. M., Weingartner J. C., Martin P. G., Volk K., Ferland G. J., 2004, in Meixner M., Kastner J. H., Balick B., Soker N., eds, Astronomical Society of the Pacific Conference Series Vol. 313, Asymmetrical Planetary Nebulae III: Winds, Structure and the Thunderbird. p. 380 (arXiv:astro-ph/0310216)

##  Appendix A Comparison with WARPFIELD-2

Figure 22: Minimum Star-formation efficiency $\epsilon_{\mathrm{SF}}$ across
the cloud mass and surface density parameter space for comparison to Fig. 5
from Rahner et al. (2019). The colormap represents an interpolation over the
data points shown.

To compare our implementation of the evolutionary model with the work carried
out by Rahner et al. (2019), we compare the minimum star formation
efficiencies needed as a function of cloud surface density and cloud mass for
a homogeneous cloud. For this comparison, we turn off the Lyman-$\alpha$
feedback. We also use the same stellar templates as the ones used in that
work. The terminal velocities for SNe ejecta are fixed to
$v_{\mathrm{sn}}=10^{4}~{}{\mathrm{m/s}}$ as mentioned in Rahner et al.
(2017), and we assume that the same is done in Rahner et al. (2019). The
parameter space used for this comparison is given in Table
LABEL:Table:benchmarking_parameter_space Our results for the minimum
$\epsilon_{\mathrm{SF}}$ required to disrupt a homogeneous cloud are shown in
Fig. 22. These results are in very good agreement with those in Rahner et al.
(2019). A more detailed comparison was not possible as WARPFIELD-2 is not
publicly available.

Table 3: The $n_{\mathrm{cl}},~{}M_{\mathrm{cl}},~{}Z,~{}\epsilon_{\mathrm{SF}}$ parameter space for the comparison with WARPFIELD-2. Parameter | Min | Max | Step size
---|---|---|---
$\log n_{\mathrm{cl}}~{}(\mathrm{~{}cm}^{-3})$ | $0.9$ | $2.7$ | $0.15$
$\log M_{\mathrm{cl}}/M_{\odot}$ | $5$ | 7 | $0.2$
$Z$ | $0.014$ | $-$ | $-$
$\epsilon_{\mathrm{SF}}~{}\%$ | 1, 1.5, 2.5, 3.5, 4.5, 6.0, 7.5

##  Appendix B Determining molecular fraction in the shell

The model described in Krumholz (2013) calculates $f_{\mathrm{mol}}$ as:

$\displaystyle f_{\mathrm{mol}}=\begin{cases}1-3s/(4+s)&\text{ for }s<2\\\
0&\text{ for }s\geqslant 2\end{cases},$ (41) $\displaystyle
s\equiv\frac{\ln\left(1+0.6\chi+0.01\chi^{2}\right)}{0.6\tau_{\mathrm{c}}},$
$\displaystyle\chi\equiv
7.2U_{\mathrm{MW}}\left(\frac{n_{\mathrm{CNM}}}{10\mathrm{~{}cm}^{-3}}\right)^{-1},$

where $U_{\mathrm{MW}}$ is the UV radiation field relative to the solar
neighborhood’s average interstellar radiation, with a flux of $3.43\times
10^{-8}$ photons s-1 cm-2 Hz-1 at $1000$ Å (Draine, 1978), and
$n_{\mathrm{CNM}}$ is the cold neutral medium’s number density.
$\tau_{\rm{c}}$ represents the cloud’s optical depth. We adapt this model to
determine the shell’s depth where the gas becomes molecular, using the stellar
SED around $1000$ Å to compute the flux at a specific shell location and
employing $n_{\rm{sh}}$ and $\tau_{\rm{d}}$ from Eqns. (16), (17) as
$n_{\mathrm{CNM}}$ and $\tau_{\rm{c}}$, respectively.

##  Appendix C Atomic Hydrogen column density trends

Figure 23: Atomic Hydrogen column density as a function of the density at the
inner face of the shell and its radius calculated using the equations in Sec.
2.1.3 for $Z=0.02$. The top row shows the results when the shells are
irradiated by a stellar cluster of age 1 Myr, while the bottom row uses a
cluster of age 8 Myr. The columns vary the mass of the shell as listed in the
columns of the top row. The mass of the stellar cluster irradiating the shells
is $5\,\times\,10^{4}\,\rm{M_{\odot}}$ in all cases. Regions that appear empty
correspond to scenarios with low inner-edge density and large radii, which are
conditions under which the shells are optically thin to ionizing radiation.
The impact of the presence of molecular Hydrogen can be seen in the reduction
in the atomic columns for large radii (lower UV flux) and high density (lower
dissociation of H2) cases, see Appendix B. Figure 24: Same as Fig. 23, but for
$Z=0.004$.

To understand the influence of various parameters on atomic Hydrogen column
densities, we apply the equations detailed in Sec. 2.1.3 while treating the
shell’s inner edge density, radius, and mass as independent parameters. This
method provides a generalized version of our models. Additionally, to examine
the influence of cluster age, we analyze a young system (1 Myr) and a system
where the most massive stars have died (8 Myr), keeping the stellar mass
irradiating the system constant at $5\,\%$ of the largest shell considered,
i.e., $M_{\star}=5\,\times\,10^{4}\,\rm{M_{\odot}}$.

Fig. 23 shows the atomic Hydrogen column density as a function of the density
at the inner face of the shell and its radius calculated for $Z=0.02$, while
Fig. 24 shows the same for $Z=0.004$.

It is clear that higher column densities of neutral Hydrogen are associated
with clouds which promote lower shell radii, higher masses, and higher inner
edge densities. In our models, while these parameters are interdependent,
specific inferences can still be made for certain types of clouds. For
instance, during the early phases of expansion, clouds with higher density
tend to have shells with smaller radii, higher swept-up mass, and higher inner
edge density compared to their low-density counterparts. This leads to higher
atomic Hydrogen column densities (if any, depending on the swept-up mass and
the ionizing radiation strength) at earlier stages, which in turn results in a
stronger coupling with Ly$\alpha$ radiation. Such a dynamic can significantly
impact the evolution of a system, particularly in environments with lower
metallicities. Likewise, clouds (and, consequently, shells) with higher masses
can result in high atomic Hydrogen column densities, even when the shell radii
are moderately large and the inner edge densities are relatively low.

##  Appendix D Absorption of ionizing radiation by dust

We estimate the fraction of ionizing photons that escape absorption by dust,
denoted as $f_{\rm{ion}}$ (or $1-f^{\rm{ion}}_{\rm{abs,\,dust}}$), using the
fitting function from Draine (2011). This analysis is valid for hydrostatic
shells, as is appropriate for our Cloudy models.

$f_{\text{ion
}}\left(\beta,\gamma,\tau_{d0}\right)\approx\frac{1}{1+(2/3+AB)\tau_{d0}}+\frac{AB\tau_{d0}}{1+B\tau_{d0}},$
(42)

with

$\displaystyle A$ $\displaystyle=\frac{1}{1+0.75\gamma^{0.65}\beta^{-0.44}},$
(43) $\displaystyle B$ $\displaystyle=\frac{0.5}{1+0.1(\gamma/\beta)^{1.5}},$
(44) $\displaystyle\beta$
$\displaystyle\equiv\frac{L_{\mathrm{n}}}{L_{\mathrm{i}}},$ (45)
$\displaystyle\gamma$
$\displaystyle=11.2\,T_{4}^{1.83}\left(\frac{18\mathrm{eV}}{\langle
h\nu\rangle_{\mathrm{i}}}\right)\left(\frac{\sigma_{d}}{10^{-21}\mathrm{~{}cm}^{2}}\right).$
(46)

Here, $\tau_{d0}$ denotes the dust optical depth within the ionized region,
$\langle h\nu\rangle_{\mathrm{i}}$ represents the mean energy of the ionizing
photons, and $T_{4}$ is the temperature in the H ii region in units of
$10^{4}\,$K.

We utilize the output from Cloudy to determine $\tau_{d0}$. In our
calculations, the ionizing region is delineated by zones where the electron
fraction exceeds $0.5$. We adopt an absorption cross section consistent with
the H ii regions’ dust model used in the Cloudy models as $\sigma_{d}$ here.
Additionally, the mass-weighted temperature is directly derived from these
models. Both $\langle h\nu\rangle_{\mathrm{i}}$ and $\beta$ are extracted
directly from the stellar templates.

##  Appendix E Line ratios of highly attenuated objects

Figure 25: Zone by zone emission from a model where the shell is deeply
embedded in the cloud.

In Fig. 13, we observe that highly attenuated objects inhabit a unique area on
the BPT diagram. This behavior is further understood by analyzing the outputs
from Cloudy. Specifically, Cloudy provides both intrinsic and emergent line-
luminosities. The intrinsic line luminosities factor in dust grains across all
processes that influence line production without taking into account
foreground attenuation. In contrast, emergent luminosities incorporate this
foreground attenuation. Notably, when we use intrinsic line luminosities, this
set of models reverts to its customary position on the BPT diagram.

Delving deeper, Fig. 25 offers a zone-by-zone emission analysis for an early
age model, at a stage where the shell remains deeply embedded within its birth
cloud. Both H$\alpha$ and [N ii] emission originating from the H ii region are
considerably attenuated—a fact underscored when one contrasts the intrinsic
with the emergent zone luminosities. Concurrently, H$\alpha$ possesses a
secondary non-recombination contribution from outside the H ii region. In an
environment devoid of dust, where the H ii region is clearly visible, this
contribution would go unnoticed. However, with the H ii region being
significantly attenuated, this secondary contribution becomes the primary
driver of the H$\alpha$ line luminosity emerging from the cloud. For the [N
ii] emission, however, the emergent contribution primarily stems from the H ii
region itself.

This phenomenon results in such models being positioned in the distinctive
area of the BPT diagram, as illustrated in both Fig. 12 and 13. A similar
explanation could also be given for the line ratios encountered in the case of
infalling shells. However, due to the pronounced attenuation, such models are
unlikely to be detected.

##  Appendix F Examples of the evolution of key shell properties

Figure 26: Evolution of key shell properties as a function of the parameters
of the model for $Z=0.02$: The top row shows the effect of changing the cloud
density of models with $\epsilon_{\mathrm{SF}}=5.0\%$, and
$\log(M_{\mathrm{cl}}/M_{\odot})=6.0$. The second row shows the effect of
changing the cloud mass of a system with $\epsilon_{\mathrm{SF}}=5.0\%$, and
$n_{\mathrm{cl}}=320.0~{}\mathrm{cm^{-3}}$. The last row depicts the impact of
changing the star-formation efficiency for a system with
$n_{\mathrm{cl}}=320.0~{}\mathrm{cm^{-3}}$, and
$\log(M_{\mathrm{cl}}/M_{\odot})=6.0$. The first, second, and third columns
show the shell radius, velocity, and inner-face density for the shells,
respectively. The fixed parameters of each row are listed in bold in the first
column, while the legend title gives the parameter being varied.

Figure 27: Same as Fig. 26, but for $Z=0.004$.