††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Engineering morphogenesis of cell clusters with differentiable programming

Ramya Deshpande  School of Engineering and Applied Sciences, Harvard
University, Cambridge MA 02138  Francesco Mottes  School of Engineering and
Applied Sciences, Harvard University, Cambridge MA 02138  Department of
Computational Biology, University of Lausanne, Lausanne, Switzerland  Ariana-
Dalia Vlad  Harvard College, Cambridge, MA 02138  Michael P. Brenner  School
of Engineering and Applied Sciences, Harvard University, Cambridge MA 02138
Alma dal Co  Department of Computational Biology, University of Lausanne,
Lausanne, Switzerland

###### Abstract

Understanding the rules underlying organismal development is a major unsolved
problem in biology. Each cell in a developing organism responds to signals in
its local environment by dividing, excreting, consuming, or reorganizing, yet
how these individual actions coordinate over a macroscopic number of cells to
grow complex structures with exquisite functionality is unknown. Here we use
recent advances in automatic differentiation to discover local interaction
rules and genetic networks that yield emergent, systems-level characteristics
in a model of development. We consider a growing tissue with cellular
interactions are mediated by morphogen diffusion, differential cell adhesion
and mechanical stress. Each cell has an internal genetic network that it uses
to make decisions based on its local environment. We show that one can
simultaneously learn parameters governing the cell interactions and the
genetic network for complex developmental scenarios, including the symmetry
breaking of an embryo from an initial cell, the creation of emergent chemical
gradients,homogenization of growth via mechanical stress, programmed growth
into a prespecified shape, and the ability to repair from damage. When
combined with recent experimental advances measuring spatio-temporal dynamics
and gene expression of cells in a growing tissue, the methodology outlined
here offers a promising path to unravelling the cellular basis of development.

††preprint: APS/123-QED

Morphogenesis, the emergence of distinct anatomical forms from dividing
cellular assemblies, hinges on the delicate coordination of diverse cellular
and molecular processes across multiple spatial and temporal scales.
Programming cells to achieve desired morphogenetic outcomes is extremely
challenging. Unlike systems guided by centralized control signals [1], most
developmental complexity arises from the interplay of local communications,
sensing, and information processing within individual cells, with minimal
external influence [2]. The emergent nature of collective cell phenomena,
coupled with the complexity of modeling them in adequate detail, imposes
severe limits to the effectiveness of computational approaches to inverse
design.

Recent experimental advances offer unprecedented opportunities for
morphogenetic engineering [3, 4], enabling the engineering of tailored
cellular interactions, with applications ranging from population control [5]
to programming multicellular assembly patterns [6]. Leveraging the intrinsic
self-organizing tendencies of cell clusters has also facilitated the
production of diverse organoids that mimic multiple aspects of real tissues
and organs [7, 8, 9]. The ability to quantitatively guide existing self-
organizing developmental programs or engineer novel ones holds significant
implications for both fundamental scientific understanding and biotechnology
[10, 11]. Many of these successes, however, heavily rely on manually crafted
and qualitative cellular interaction rules and protocols. Heuristic
experimental rules often exhibit drawbacks, including system specificity,
limited robustness, and insufficient reproducibility [12]. Developing more
versatile and efficient tools for biological design has the potential to
extend the reach and impact of experimental techniques.

Inverse design tasks are usually translated into optimization problems,
wherein the objective is to minimize the discrepancy between a property of the
simulation system state and its desired value. While sophisticated
computational models of cell behavior exist [13, 14], their applicability to
reverse engineering remains limited, particularly in high-dimensional
parameter spaces where parameter sweeps and evolutionary searches become
inefficient. Automatic differentiation (AD) algorithms, on the other hand,
enable efficient calculation of sensitivities of observables with respect to
input parameters. While originally formulated for training large neural
networks [15, 16], the same infrastructure can be used to differentiate
through simulations of physical processes, enabling gradient-based
optimization as a viable alternative to tackle hard inverse problems of
various scientific interest. Examples include the design of self-assembling
materials [17], the establishment of data-driven discretizations for partial
differential equations [18], nonequilibrium statistical mechanics [19] and the
design of quantum devices [20]. Neural Cellular Automata [21] represented an
initial demonstration of how these methods could be used to engineer
collective behaviors inspired by biology in lattice models. Hiscock, instead,
showed that AD could be used to successfully learn meaningful coupling in ODE
models of genetic networks [22].

Here, we consider a model of a tissue composed of cells, where each cell can
divide, grow, exert stress, as well as excrete and sense morphogens. Cellular
decisions are made in response to the local environment, via an internal
genetic network (Fig. 1). Starting from random initial parameters, we
demonstrate that we can learn the parameters governing both cellular and
genetic interactions so that the cluster of cells undergoes non-trivial
developmental scenarios, ranging from symmetry breaking of an embryo from a
progenitor cell, to the creation of emergent chemical gradients, to
homogenization of growth via chemical signalling and mechanical stress
response, to programmed growth into a pre-specified shape, which is able to
repair itself from damage. The framework allows for the probing of the
relationship between the capabilities of individual cells and the emergent
properties of the resulting tissues.

Figure 1: Modeling and Optimization Framework. (a) Model overview. Decision-
making circuits in cells sense local chemical and mechanical inputs to
regulate secretion, division probability, and adhesion. The environment
responds to cellular outputs by handling chemical diffusion, cell division,
and mechanical relaxation. (b) Input signals are integrated and elaborated in
each cell by internal gene regulatory networks. Green gene nodes sense inputs,
blue genes are intermediates and red genes can modify cellular actions. Gene
expression levels in each node are maintained from one simulation step to the
next, providing memory to the system. (c) Gradient calculation and
optimization. Following the REINFORCE approach, the log probability gradient
of each time step is weighted by action rewards to estimate the average
gradient and parameters are updated through gradient descent. Gradients flow
through the whole differentiable simulation, accounting for environment
contributions.

## Model and Optimization Framework

### Forward Model

We model tissue growth using soft spheres interacting through a Morse
potential [23]. We mostly consider two-dimensional systems since three-
dimensional systems rapidly become too computationally intensive to optimize
with our current implementation. Our methods are nonetheless applicable in
principle to any physical model of cell interactions and we show some simple
experiments with 3D systems. In our model, cells grow at a fixed rate until
reaching a maximum radius, leading to continuous mechanical reorganization.
Cells secrete $M$ morphogens, measure local concentrations and spatial
gradients of these chemicals, and sense mechanical stress. This locally sensed
information guides decisions on chemical secretion modulation and division
probability (Fig. 1a). Cell division events are stochastic, producing two
daughter cells, each with half the volume of a fully grown mother cell.
Decision-making circuits are modeled as regulatory networks comprising $N$
interacting genes with trainable excitatory or inhibitory couplings (Fig. 1b).
The regulatory network processes the information sensed by individual cells to
output rates of division, chemical secretion, etc. In most examples, the
regulatory network weights are shared by all cells, though we also explore one
scenario where two cell types have partially different regulatory networks.
Further mathematical detail can be found in the Supplementary Information.

### Inverse Design

Our objective is to discover local growth rules that enable cells to develop
into a final configuration with desired properties. To effectively navigate
the high-dimensional, non-convex parameter space of local interaction rules,
we employ gradient descent schemes, widely used in modern machine learning for
such optimizations. Our entire simulation is written to be automatically
differentiable, leveraging recent advances from the deep learning and
reinforcement learning communities. We utilize the JAX library [24] for
efficient numerical and automatic differentiation algorithms, and other JAX-
based libraries such as Equinox [25] and JAX-MD [26] for molecular dynamics.

We employ score-based methods such as REINFORCE [27] to address the intrinsic
stochasticity of growth dynamics, which implies that the loss function is not
necessarily differentiable (Fig. 1c). In the REINFORCE update, rewards and
penalties are assigned to each action taken (i.e. each cell division) to
provide the weighting signal for the gradient of the trajectory probability.
This weighted gradient is then used to nudge the weights of the gene network
in the right direction. In practice, the system gradually learns which
division events are the most favorable and increases their probability for the
next simulation. To guide the reader possibly unfamiliar with these concepts,
we also include a set of in-depth tutorial notebooks in the GitHub repository
[28] and a longer theoretical explanation in the Supplementary Information.
Since the dynamics of the environment, including chemical secretion and
diffusion, depend partly on the same parameters (the gene network couplings)
that cells use to decide whether to divide or not (their “policy”), it is
essential to evaluate the gradient of the environment dynamics with respect to
the gene network parameters. This necessity makes automatic differentiation
particularly crucial in this context. Once the necessary gradients are
calculated, optimization is performed using well-established gradient descent
schemes, such as Adam [29]. The learned networks are then further simplified
by removing all edges with weights smaller than an experiment-dependent
threshold, such that the functional backbone of pruned network be highlighted
while preserving the optimized properties.

Figure 2: Chemical Regulation of Cellular Adhesion. Learned chemical
regulation of cellular adhesion guides the assembly of different spatial
cluster arrangements. The top part of the panels show the training loss curves
for three different optimization targets. The insets show the final state of a
simulation after the learning, with colors denoting different cell types. The
bottom sections show the learned and pruned gene network; Cad1, Cad2, and Cad3
are the homotypic cadherins of the three cell types, and Het is the
heterotypic cadherin expressed in all cell types. (a) Lattice-like structure.
Cells are arranged to minimize the sum of pairwise distances between
heterotypic cells. (b) Lobed structure. Cells minimize the sum of pairwise
distances between homotypic cells. (c) Core-shell structure. Cells are
arranged to minimize the distance of cells from a predefined radius, dependent
on their type. The final configuration is here shown as a cross-section to
highlight the internal structure.

## Chemical Regulation of Cellular Adhesion

The differential adhesion hypothesis posits that populations of cells can
adopt specific morphologies based on differences in adhesive strength [30].
During morphogenesis, a cadherin code regulates the spatial organization of
developing tissue [31, 32]. Engineered cellular circuits can use gene
regulatory networks to dynamically modulate adhesion and program self-
organization [6]. Here, we show that gene interactions can be learned that use
chemical signaling to modulate homotypic and heterotypic adhesion to achieve a
target spatial pattern in a population with different cell types.

In our model, a population of three cell types is subject to Brownian motion
and interacts via a Morse pairwise potential. The well depth of the potential
between a pair of cells is dictated by the concentration of homotypic and
heterotypic “cadherins” in each cell. Each cell type produces its respective
homotypic cadherin and a shared heterotypic cadherin, regulated by chemical
signaling. To achieve a desired spatial organization, we constructed three
loss functions based on pairwise distances between cells. For a lattice-like
structure, we minimized the sum of pairwise distances of cells of different
types (Fig 2a). For organizing the structure into three lobes, we minimized
the sum of pairwise distances of cells of the same type (Fig 2b). To achieve a
core-shell structure, we minimized the distance from each cell to a pre-
specified radius based on the cell type (Fig 2c).

In each case, the optimization learned simple mechanisms of spatial control.
The optimization process identified effective solutions by allowing cells to
observe their own type and regulate the pairwise potential between cell types.
For a lattice-like structure, cells upregulated their heterotypic cadherin and
downregulated their homotypic cadherin. Conversely, to create a lobed
structure, they upregulated their homotypic cadherin and downregulated their
heterotypic cadherin. These patterns of homo- and heterotypic adhesion are
consistent with observations from cell culture sorting experiments [31].
Achieving core-shell structures instead required a more complex regulatory
network (Fig 2c).

## Spatial Control of Tissue Growth

Understanding how to program the growth of complex shapes can aid in our
comprehension of developmental body plans driven by chemical patterning.
Rationally designed synthetic gene circuits have already been shown to be able
to drive the formation of target multicellular structures found in
morphogenesis [6, 33]. We investigate here how to use our framework to
discover new decision-making circuits that utilize local chemical interactions
to grow pre-specified shapes. These rules are, as above, formulated as gene
regulatory circuits. We focus in particular on two fundamental shapes commonly
found in developmental biology: elongation and branching. The symmetry-
breaking elongation process is the first geometric-defining process to happen
in many developmental programs and is vital, for example, for establishing the
antero-posterior axis in animals [34]. Branching is possibly the next most
ubiquitous process in the development of multicellular organisms and branched
structures are found from plant roots to blood vessels and mammary epithelial
cultures [35].

We show below that gene regulations that achieve these final structures can be
learned in our model. Notably, we also show that mechanisms optimized solely
to achieve a specific final structure exhibit emergent properties that they
were not explicitly trained to achieve. In particular, gene circuits learn to
leverage, modify and generate chemical gradients and they display basic
regenerative capabilities.

Figure 3: Elongation. Maximization of the squared sum of the z-coordinates of
cells, resulting in vertical elongation in the optimized cluster. (a) Loss
curve for one successful optimization run. On the left the final configuration
of a simulation with randomly initialized parameters, on the right the final
state after learning. (b) Visualization of the optimized and pruned gene
network, displaying activation and inhibition interactions that regulate cell
behavior during elongation. (c) Two examples of emerging chemical gradients
that are used to guide growth.

### Elongation

We consider a three dimensional model of dividing cells that we optimize to
minimize the squared sum of the z-coordinates of cells, thus requiring the
cluster to elongate vertically (Fig. 3a). Local chemical concentrations, local
estimated chemical gradients, cell size, and division propensities are passed
as inputs into the cells’ gene network (Fig. 3b), that infer positional
information through local physical communication. Simulations start from a
single cell and grow to a total of 120 cells.

In order to globally assume an elongated shape, cells learn to create emergent
chemical gradients or control existing ones, even though this task is not
included in the optimization objective (Fig 3c). The cells learn to divide
primarily at one or both ends of the growing cluster. Since the cluster edge
itself moves as a consequence of new cells being added, the gene networks must
implement a non-trivial control mechanism to keep the cluster growing.

Figure 4: Branching. Optimizing growth in a V-shaped reward landscape results
in a learned branched structure of the final cluster. (a) Loss curve for one
successful optimization run. On the left, the final configuration of a
simulation with randomly initialized parameters. On the right, the final state
after learning with a visualization of the reward mask (red areas are
associated to penalties, white ones with rewards). (b) Two examples of
emerging chemical gradients that are used to guide growth. (c) Final optimized
state after 129 cell divisions, the same number used for the training (top)
and with 500 extra cell divisions. Note the spontaneous branching on the main
arms. (d) Loss of the final state after removal, at random, of a fraction of
cells. The red dashed line is the loss associated to the final simulation
state before cell removals.

### Branching

We simulate the growth of a cluster from a single cell to 130 cells,
restricting ourselves to the 2-dimensional case for computational simplicity.
With a similar spirit to diffusion-based geometry [35], we simulate permeable
diffusive boundaries between the cell cluster and the surrounding implicit
medium by connecting boundary cells to a ghost chemical sink node. Gene
networks have access to the local chemical concentrations, estimates of local
chemical gradients, their own cell radius and division propensity. They do not
have access to their own absolute spatial coordinates.

To learn more complex branched structures, we construct a position-based shape
loss. At each time point, cell divisions that place cells within the region
designating the desired shape are assigned a reward, the others a penalty. In
this case, we impose a V-shaped spatial reward structure (Fig. 4a on the
right) augmented with an asymmetry penalty on the x-coordinate, in order to
promote the growth of both arms at the same time. Also in this case, a variety
of chemical patterns emerge in the learned simulation (Fig 4b), and the
cluster learns to divide cells primarily at the tips of the branching
structure. Clusters visually acquire a branched shape, even when simulated for
500 cell divisions more than optimized for. While shape errors clearly emerge,
they do mostly in the form of smaller branches that begin to sprout from the
main branched structure (Fig 4c). Furthermore, the learned network displays
emergent regenerative capabilities. After reaching the final branched
configuration, we remove a fraction of cells chosen at random and we allow the
cluster grow back to its original size. As shown in Fig 4d, eliminating up to
60% of cells still results in the regeneration of a branched structure. Even
more interestingly, it appears that randomly removing cell actually promotes
error-correction in the regrowth, reaching even slightly better performances
in the training loss.

Figure 5: Mechanical Mechanisms for Homogeneous Growth. (a) Schematic of setup
of optimization problem for learning uniform growth and definition of loss
function. (b) Visualization of the final configuration of the cluster from an
optimized and unoptimized simulation. The unoptimized simulation does not
regulate division propensities based on local inputs and thus, growth reflects
the growth factor gradient. Cells are colored by the division propensity. (c)
Visualization of trained gene network. (d) Gene network learned with a power
law morphogen field readily works on three other kinds of morphogen fields: an
exponential field, a field with two peaks instead of one, and a field with an
off center peak. Losses are shown for each field type, averaged over 100
simulations. (e) Division rate divided by the morphogen field concentration is
plotted for each kind of simulation from (d). This shows how cells regulate
division probabilities without the influence of the field; results show that
the way cells regulate division adapts to the kind of field due to the learned
stress-based mechanism.

## Mechanical Control of Cell Proliferation

Mechanical signaling has long been recognized as an important component of
morphogenesis. In individual cells, mechanosensation influences gene
expression through several identified biochemical pathways that can transduce
forces [36]. On a tissue-level scale, mechanical signals can communicate local
instructions to build a global shape [36, 37, 38, 39]. In tissues, mechanical
interactions can act on larger length scales than biochemical signals due to
the propagation of tension through connected cells [39]. The buildup of local
mechanical stress, additionally, can influence global growth of a developing
structure. As an example, mechanical feedback has been posited as a mechanism
for growth arrest in the wing imaginal disc of Drosophila [40, 41]. Here, we
explore these ideas by optimizing over a model inspired by Hufnagel et al. to
learn mechanisms that couple mechanics to growth.

In our model, a radial growth factor gradient scales a division propensity
calculated by gene networks in each cell, such that cells closer to the origin
have a higher probability of division than cells at the edges (Fig. 5a). The
growth factor gradient follows a power law and is fixed with respect to the
simulation. We run the simulation for 200 cell division steps starting from
one cell, with each cell modulating its secretion of two morphogens and its
division propensity. We seek local growth rules that achieve homogeneous
growth (homogeneous division propensities) throughout the cluster despite the
nonuniform effect of the growth factor. The gene network that calculates the
propensity of division does not see the growth factor as an input; therefore,
cells have to infer this information by sensing stress or communicating
through morphogens. The learned gene network can sense morphogen
concentrations, morphogen gradients, and stress as cellular inputs. We
optimize couplings between genes that sense these inputs, four internal genes,
and three output genes that modulate secretion rates of the two morphogens as
well as the division propensity. The optimization minimizes the coefficient of
variation of division propensities at the end of the simulation:

$\mathcal{L}=\text{CV}[\textbf{d}(T)]=\frac{\sigma[\textbf{d}(T)]}{\mu[\textbf{d}(T)]}$

where $\textbf{d}(T)$ are the division propensities at final time $T$.

Simulations run using the optimized model show a clear separation between the
final costs of simulations run using untrained versus trained parameters (Supp
Fig. S3b). Visualization of simulations confirms that trained parameters
resulted in homogeneous growth (Fig. 5b). The learned gene circuit utilized
both chemical and mechanical information to modulate division propensity (Fig.
5c, Supp Fig. S3d). To probe the robustness of the learned mechanism, we apply
the trained gene circuit to simulations with other types of growth factor
gradients: (a) an exponentially decaying growth factor gradient, (b) a
gradient with two symmetric peaks, and (c) a gradient with a single asymmetric
peak. Although we did not show any of these morphogen patterns during
training, the gene circuit optimized on the power law gradient is still able
to homogenize growth for the other types of gradients (Fig. 5d). Furthermore,
the outputs of the gene network seem to track the shape of the fields, as well
as patterns of stress accumulation (Fig. 5e, Supp Fig. S3c), consistently with
a mechanical feedback mechanism. This is reminescent of the mechanism proposed
by Hufnagel et al., where nonuniform growth leads to a buildup of mechanical
stress at the center of the cluster, providing a feedback signal that cells
can use to regulate their growth rates.

## Chemical Regulation of Tissue Homeostasis

Maintaining stable ratios of various cell populations in tissue is important
for proper functioning, and disruption of this can lead to disease [42]. To
stably maintain homeostasis, the underlying cellular circuitry must prevent
runaway growth for either cell type. Cells communicate via growth factors, and
cell circuits that receive and send growth factors maintain tissue
homeostasis, sometimes regenerating after perturbations. For example, co-
cultures of macrophages and fibroblasts can reach stable ratios, even if
initially plated with ratios that differ 100-fold [43]. Here, we aim at
learning chemical mechanisms that can maintain a 50/50 ratio of two different
cell types in a growing cell cluster.

In our model, two cell types are instantiated, each secreting one of two
different chemicals. The simulation is seeded with an imbalanced initial state
of $N_{init}=20$ total cells and $n_{r}=4$ red cells, and uniform division
propensities. After 100 cell divisions, this results in a final cluster that
preserves the initial imbalance (Fig 6a left). Gene interaction couplings are
then optimized to minimize the number difference between the two cell types
during growth (Fig 6a right). The loss at time $t$ takes the form:

$\mathcal{L}_{t}=|n_{r}(t)-n_{b}(t)|$

where $n_{r}$ and $n_{b}$ are the number of red and blue cells respectively,
with $n_{r}+n_{b}=N$ being the total number of cells at time $t$.

Cells can take division decisions only based on local communication of
chemicals, and therefore they must infer information about population
differences from chemical interactions. The optimization procedure
successfully learns mechanisms that decrease the initial imbalance – such as
the one depicted in Fig. 6b – resulting in equal proportions of cell types.
Fig. 6c shows instead the chemical landscape associated with the learned final
state, where one can observe the chemicals being associated with the two
different cell types. Remarkably, the learned gene regulation displays
robustness and generalization properties that were not explicitly optimized
for. Fixing the initial imbalance to the one used during training
($n_{r}/N=.2$), even though the training simulations only go through 100 cell
divisions, the learned regulation maintains homeostasis also for simulations
that run for much longer with no sign of efficiency degradation (Fig. 6d). If
we instead fix the number of cell divisions to 100 and change the initial
state imbalance, we observe in Fig. 6e that homeostasis is reasonably
maintained for up to almost twice the initial imbalance. After this threshold
the circuit starts to over-compensate, and loses its advantage over completely
random cell divisions.

Figure 6: Chemical Regulation of Homeostasis. (a) An initially imbalanced
state (in the middle) grows into an equally unbalanced final state with
untrained network parameters (left) and into a balanced state with learned
parameters (right) (b) Visualization of the learned gene network. (c)
Visualization of the chemical concentrations at the final state of the
optimized simulation. (d) Average over 100 realizations of final proportion of
red cells at the end of simulations of increasing lengths. (e) Average over
100 realizations of final proportion of red cells at the end of simulations
with increasing proportions of red cells in the starting state. In both plots
the shaded area corresponds to average plus or minus one standard deviation.

## Discussion

Our results demonstrate the power of gradient-based optimization methods in
discovering mechanisms within the space of physical, cell-based models for
collective behavior. By taking gradients through a physical model of cells
interacting written using differentiable code, we were able to learn decision-
making programs within individual cells to orchestrate growth into target
shapes or with target properties. Additionally, we were able to overcome
challenges associated with differentiating through the discrete nature of cell
division events and non-differentiable loss functions through reinforcement
learning based approaches. The mechanisms that were learned took the form of
gene regulatory networks, which readily offers a biological interpretation
that could even inform experiments. The learned mechanisms were local in
nature; cells did not have access to global positional information and had to
instead infer it by sensing chemicals and mechanics and communicating with
nearby cells through these physical channels.
Our results are a promising first step towards inferring mechanisms in
existing multicellular systems or learning how to program them for desired
behaviors. In particular, morphogenetic mechanisms frequently utilize local
feedback to coordinate intricate forms, often through chemical and/or
mechanical signals. Examples include a myriad of developmental systems, from
the development of the wing imaginal disc in Drosophila to the periodic
segmentation of blocks in somitogenesis. The development of robust and useful
organoid systems also requires an understanding of spatiotemporal chemical and
physical cues needed to guide the growth of stem cells. Physical models could
be optimized to reproduce experimental data from time-lapse imaging of such
developmental systems to gain insights into individual decision making
functions that recapitulate complex behaviors. Since we have shown that gene
regulatory networks can be programmed to drive complex self-organization,
further incorporation of experimental knowledge or distillation of the learned
networks could even be instantiated experimentally in synthetic circuits.
Furthermore, the soft sphere model we have presented here is simply a choice;
the physical details can be modified according to the biological context the
model is designed to ask questions for. The optimization can also be adjusted
accordingly depending on the question being asked: mechanical parameters could
be optimized over to learn regulation of cellular potentials, or
hyperparameters like the number of chemicals. Questions about minimal
mechanisms, robustness to variation and perturbations, or growth principles
could also be framed and investigated. Gene expression datasets could be
incorporated to bridge the gap between known and unknown information in
biological systems. Our work opens up a multitude of avenues for further
explorations of the landscape of growth programs for desired emergent
behaviors in living systems.

## Materials and Methods

Fully differentiable, hardware-accelerated simulations were coded using the
JAX python library [24] and other libraries in the JAX ecosystem. In
particular, JAX-MD [26] was used to perform molecular dynamics simulations for
cell mechanics and Equinox [25] as the simulation backbone.

We proceed to give a very succinct overview of the main components of our
model. A full mathematical account of simulation and optimization details is
given in the Supplementary Information.

### Cell division and growth

Each cell $i$ is endowed with a division propensity $d_{i}$ calculated by the
internal gene network. One cell divides at each time step, with probability:

$P(i\ \text{divides})=\frac{d_{i}}{\sum_{i}d_{i}}$

The division substitutes one mature cell with two newborn cells with half the
volume. The angular orientation of the division plan is chosen uniformly at
random.

Cells’ radii $R_{i}$ grow linearly in time up to a maximum radius $R_{max}$:

$R_{i}(t+\Delta t)=\min(R_{i}(t)+\Delta R,R_{\text{max}})$

### Mechanical interactions

Cells are modeled as adhesive soft spheres. Interactions are defined by the
Morse energy potential:

$V_{ij}(r)=\epsilon_{ij}\left[1-e^{-\alpha_{ij}(r_{ij}-\sigma_{ij})}\right]^{2}-\epsilon_{ij}$

where $\epsilon_{ij}$ is the well depth for the potential between particle $i$
and $j$, $\alpha_{ij}$ is the range of interaction, $\sigma_{ij}$ is the sum
of radii of particles $i$ and $j$, and $r_{ij}$ is the distance between the
two particles.

### Diffusion

We assume that chemicals diffuse freely in the environment and are taken up by
cells at a constant rate. It is also assumed that diffusion is much faster
than all other processes and can be therefore be considered at steady state.
For any one chemical $k$:

$\frac{\partial c_{k}}{\partial t}=D_{k}\nabla^{2}c_{k}-K_{k}c_{k}+S_{k}=0$

The equation is discretized and solved on the lattice formed by the cells’
centers with the graph Laplacian technique.

### Genetic regulatory interactions

We model internal cell decision functions with a simple ODE model of gene
regulation inspired by Hiscock [22]:

$\frac{dg_{i}}{dt}=\phi(\Sigma_{j}W_{ij}g_{j})+I_{i}-k_{i}g_{i}$

Considering the $i$-th gene, $g_{i}$ is the expression level, $W_{ij}$ is the
regulatory influence of gene $j$ on $i$, $k_{i}$ is the degradation rate, and
$I_{i}$ is a forcing signal coming from one of the sensed cell inputs. $\phi$
is a sigmoid function that ensures positive transcription rates saturating at
high input levels.

### Gradient calculation

The gradient of simulations in Fig. 2 is calculated with standard automatic
differentiation. The gradients of all the other simulations involve dealing
with cell division events. Following REINFORCE [27], we can estimate the
gradient of the expected loss for a single simulation:

$\nabla_{\theta}\mathcal{L}(\theta)\approx\sum_{t=0}^{T-1}L_{t}\nabla_{\theta}\log\pi_{\theta}(a_{t}|s_{t})$

where $L_{t}$ is the discounted future loss, $\theta$ represents the trainable
gene network parameters and $\pi$ is the probability of the division event
parametrized by $\theta$.

###### Acknowledgements.

This work is dedicated to the memory of Alma Dal Co. We thank for insightful
discussions: the members of the Brenner Group, the members of the former Dal
Co group, Boris Shraiman and Dillon Cislo. R.D. thanks Hervé Turlier for
helpful discussions. This work was supported by the Office of Naval Research
through grant numbers ONR N00014-17-1-3029, ONR N00014-23-1-2654 and the NSF
AI Institute of Dynamic Systems (2112085).

## Data Availability

The custom library written to conduct this study is available at
https://github.com/fmottes/jax-morph/, along with tutorial notebooks.

## References

  * Rubenstein  _et al._ [2014] M. Rubenstein, A. Cornejo, and R. Nagpal, Programmable self-assembly in a thousand-robot swarm, Science 345, 795 (2014).
  * Collinet and Lecuit [2021] C. Collinet and T. Lecuit, Programmed and self-organized flow of information during morphogenesis, Nature Reviews Molecular Cell Biology 22, 245 (2021).
  * Davies [2017] J. Davies, Using synthetic biology to explore principles of development, Development (Cambridge, England) 144, 1146 (2017).
  * Velazquez  _et al._ [2018] J. J. Velazquez, E. Su, P. Cahan, and M. R. Ebrahimkhani, Programming Morphogenesis through Systems and Synthetic Biology, Trends in biotechnology 36, 415 (2018).
  * Ma  _et al._ [2022] Y. Ma, M. W. Budde, M. N. Mayalu, J. Zhu, A. C. Lu, R. M. Murray, and M. B. Elowitz, Synthetic mammalian signaling circuits for robust cell population control, Cell 185, 967 (2022).
  * Toda  _et al._ [2020] S. Toda, W. L. McKeithan, T. J. Hakkinen, P. Lopez, O. D. Klein, and W. A. Lim, Engineering synthetic morphogen systems that can program multicellular patterning, Science (New York, N.Y.) 370, 327 (2020).
  * Lancaster  _et al._ [2017] M. A. Lancaster, N. S. Corsini, S. Wolfinger, E. H. Gustafson, A. Phillips, T. R. Burkard, T. Otani, F. J. Livesey, and J. A. Knoblich, Guided self-organization and cortical plate formation in human brain organoids, Nature biotechnology 35, 659 (2017).
  * Lewis  _et al._ [2021] A. Lewis, R. Keshara, Y. H. Kim, and A. Grapin-Botton, Self-organization of organoids from endoderm-derived cells, Journal of Molecular Medicine 99, 449 (2021).
  * Brassard and Lutolf [2019] J. A. Brassard and M. P. Lutolf, Engineering Stem Cell Self-organization to Build Better Organoids, Cell Stem Cell 24, 860 (2019).
  * Sthijns  _et al._ [2019] M. M. Sthijns, V. L. LaPointe, and C. A. van Blitterswijk, Building Complex Life Through Self-Organization, Tissue Engineering. Part A 25, 1341 (2019).
  * Werner  _et al._ [2017] S. Werner, H. T.-K. Vu, and J. C. Rink, Self-organization in development, regeneration and organoids, Current Opinion in Cell Biology Cell Architecture, 44, 102 (2017).
  * Hofer and Lutolf [2021] M. Hofer and M. P. Lutolf, Engineering organoids, Nature Reviews Materials 6, 402 (2021).
  * Osborne  _et al._ [2017] J. M. Osborne, A. G. Fletcher, J. M. Pitt-Francis, P. K. Maini, and D. J. Gavaghan, Comparing individual-based approaches to modelling the self-organization of multicellular tissues, PLoS Computational Biology 13, e1005387 (2017).
  * Runser  _et al._ [2024] S. Runser, R. Vetter, and D. Iber, Simucell3d: three-dimensional simulation of tissue mechanics with cell polarization, Nature Computational Science 4, 299 (2024).
  * Anil  _et al._ [2023] R. Anil, A. M. Dai, O. Firat, M. Johnson, D. Lepikhin, A. Passos, S. Shakeri, E. Taropa, P. Bailey, Z. Chen, _et al._ , Palm 2 technical report, arXiv preprint arXiv:2305.10403 (2023).
  * Brown  _et al._ [2020] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, _et al._ , Language models are few-shot learners, Advances in neural information processing systems 33, 1877 (2020).
  * Goodrich  _et al._ [2021] C. P. Goodrich, E. M. King, S. S. Schoenholz, E. D. Cubuk, and M. P. Brenner, Designing self-assembling kinetics with differentiable statistical physics models, Proceedings of the National Academy of Sciences 118, e2024083118 (2021).
  * Bar-Sinai  _et al._ [2019] Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner, Learning data-driven discretizations for partial differential equations, Proceedings of the National Academy of Sciences 116, 15344 (2019).
  * Engel  _et al._ [2023] M. C. Engel, J. A. Smith, and M. P. Brenner, Optimal control of nonequilibrium systems through automatic differentiation, Physical Review X 13, 041032 (2023).
  * Vargas-Hernández  _et al._ [2021] R. A. Vargas-Hernández, R. T. Q. Chen, K. A. Jung, and P. Brumer, Fully differentiable optimization protocols for non-equilibrium steady states, New Journal of Physics 23, 123006 (2021).
  * Mordvintsev  _et al._ [2020] A. Mordvintsev, E. Randazzo, E. Niklasson, and M. Levin, Growing Neural Cellular Automata, Distill 5, e23 (2020).
  * Hiscock [2019] T. W. Hiscock, Adapting machine-learning algorithms to design gene circuits, BMC Bioinformatics 20, 214 (2019).
  * Koyama  _et al._ [2023] H. Koyama, H. Okumura, A. M. Ito, K. Nakamura, T. Otani, K. Kato, and T. Fujimori, Effective mechanical potential of cell–cell interaction explains three-dimensional morphologies during early embryogenesis, PLOS Computational Biology 19, e1011306 (2023), publisher: Public Library of Science.
  * Bradbury  _et al._ [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018).
  * Kidger and Garcia [2021] P. Kidger and C. Garcia, Equinox: neural networks in JAX via callable PyTrees and filtered transformations, Differentiable Programming workshop at Neural Information Processing Systems 2021 (2021).
  * Schoenholz and Cubuk [2020] S. S. Schoenholz and E. D. Cubuk, Jax m.d. a framework for differentiable physics, in  _Advances in Neural Information Processing Systems_, Vol. 33 (Curran Associates, Inc., 2020).
  * Sutton and Barto [2018] R. S. Sutton and A. G. Barto,  _Reinforcement Learning: An Introduction_, 2nd ed. (The MIT Press, 2018).
  * jax [2024] https://github.com/fmottes/jax-morph/tree/eqx/tutorials (2024), accessed: 2024-07-06.
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization (2017), arXiv:1412.6980 [cs].
  * Steinberg [1970] M. S. Steinberg, Does differential adhesion govern self-assembly processes in histogenesis? Equilibrium configurations and the emergence of a hierarchy among populations of embryonic cells, Journal of Experimental Zoology 173, 395 (1970), _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jez.1401730406.
  * Tsai  _et al._ [2020] T. Y.-C. Tsai, M. Sikora, P. Xia, T. Colak-Champollion, H. Knaut, C.-P. Heisenberg, and S. G. Megason, An adhesion code ensures robust pattern formation during tissue morphogenesis, Science 370, 113 (2020), publisher: American Association for the Advancement of Science.
  * Bao  _et al._ [2022] M. Bao, J. Cornwall-Scoones, E. Sanchez-Vasquez, A. L. Cox, D.-Y. Chen, J. De Jonghe, S. Shadkhoo, F. Hollfelder, M. Thomson, D. M. Glover, and M. Zernicka-Goetz, Stem cell-derived synthetic embryos self-assemble by exploiting cadherin codes and cortical tension, Nature Cell Biology 24, 1341 (2022), publisher: Nature Publishing Group.
  * Guo  _et al._ [2021] Y. Guo, M. Nitzan, and M. P. Brenner, Programming cell growth into different cluster shapes using diffusible signals, PLOS Computational Biology 17, e1009576 (2021).
  * Darras  _et al._ [2018] S. Darras, J. H. Fritzenwanker, K. R. Uhlinger, E. Farrelly, A. M. Pani, I. A. Hurley, R. P. Norris, M. Osovitz, M. Terasaki, M. Wu, J. Aronowicz, M. Kirschner, J. C. Gerhart, and C. J. Lowe, Anteroposterior axis patterning by early canonical wnt signaling during hemichordate development, PLoS Biology 16, e2003698 (2018).
  * Iber and Menshykau [2013] D. Iber and D. Menshykau, The control of branching morphogenesis, Open Biology 3, 130088 (2013).
  * Agarwal and Zaidel-Bar [2021] P. Agarwal and R. Zaidel-Bar, Mechanosensing in embryogenesis, Current Opinion in Cell Biology Cell Architecture, 68, 1 (2021).
  * Irvine and Shraiman [2017] K. D. Irvine and B. I. Shraiman, Mechanical control of growth: ideas, facts and challenges, Development 144, 4238 (2017).
  * Trinh  _et al._ [2021] D.-C. Trinh, J. Alonso-Serra, M. Asaoka, L. Colin, M. Cortes, A. Malivert, S. Takatani, F. Zhao, J. Traas, C. Trehin, and O. Hamant, How Mechanical Forces Shape Plant Organs, Current biology: CB 31, R143 (2021).
  * LeGoff  _et al._ [2013] L. LeGoff, H. Rouault, and T. Lecuit, A global pattern of mechanical stress polarizes cell divisions and cell shape in the growing Drosophila wing disc, Development 140, 4051 (2013).
  * Shraiman [2005] B. I. Shraiman, Mechanical feedback as a possible regulator of tissue growth, Proceedings of the National Academy of Sciences 102, 3318 (2005), publisher: Proceedings of the National Academy of Sciences.
  * Hufnagel  _et al._ [2007] L. Hufnagel, A. A. Teleman, H. Rouault, S. M. Cohen, and B. I. Shraiman, On the mechanism of wing size determination in fly development, Proceedings of the National Academy of Sciences 104, 3835 (2007), publisher: Proceedings of the National Academy of Sciences.
  * Hicks-Berthet  _et al._ [2021] J. Hicks-Berthet, B. Ning, A. Federico, A. Tilston-Lunel, A. Matschulat, X. Ai, M. E. Lenburg, J. Beane, S. Monti, and X. Varelas, Yap/taz inhibit goblet cell fate to maintain lung epithelial homeostasis, Cell Reports 36, 109347 (2021).
  * Zhou  _et al._ [2018] X. Zhou, R. A. Franklin, M. Adler, J. B. Jacox, W. Bailis, J. A. Shyer, R. A. Flavell, A. Mayo, U. Alon, and R. Medzhitov, Circuit Design Features of a Stable Two-Cell System, Cell 172, 744 (2018), publisher: Elsevier.