# AI Research Associate for Early-Stage Scientific Discovery

Morad Behandish1, John T. Maxwell III1, Johan de Kleer1
Corresponding Author, e-mail<EMAIL_ADDRESS>
###### Abstract

Artificial intelligence (AI) has been increasingly applied in scientific
activities for decades; however, it is still far from an insightful and
trustworthy collaborator in the scientific process. Most existing AI methods
are either too simplistic to be useful in real problems faced by scientists or
too domain-specialized (even dogmatized), stifling transformative discoveries
or paradigm shifts. We present an AI research associate for early-stage
scientific discovery based on (a) a novel minimally-biased ontology for
physics-based modeling that is context-aware, interpretable, and generalizable
across classical and relativistic physics; (b) automatic search for viable and
parsimonious hypotheses, represented at a high-level (via domain-agnostic
constructs) with built-in invariants, e.g., postulated forms of conservation
principles implied by a presupposed spacetime topology; and (c) automatic
compilation of the enumerated hypotheses to domain-specific, interpretable,
and trainable/testable tensor-based computation graphs to learn
phenomenological relations, e.g., constitutive or material laws, from sparse
(and possibly noisy) data sets.

## Introduction

Data-driven AI methods have been applied extensively in the past few decades
to distill nontrivial physics-based insights (scientific discovery) and to
predict complex dynamical behavior (scientific simulation) (Stevens et al.
2020). Notwithstanding their effectiveness and efficiency in classification,
regression, and forecasting tasks, statistical learning methods can hardly
ever evaluate the soundness of a function fit, explain the reasons behind
observed correlations, or provide sufficiently strong guarantees to replace
parsimonious and explainable scientific expressions such as differential
equations (DE). Hybrid methods such as constructing “physics-
informed/inspired/guided” architectures for neural nets and loss functions
that penalize both predication and DE residual errors (Raissi, Perdikaris, and
Karniadakis 2019; Wei and Chen 2019; Daw et al. 2020) and graph-nets based on
control theory and combinatorial structures (Cranmer et al. 2019; Seo and Liu
2019; Sanchez-Gonzalez et al. 2020) are all important steps towards
explainability; however, the built-in ontological biases in most machine
learning (ML) frameworks prevent them from thinking outside the box to
discover not only the known-unknowns, but also unknown-unknowns, during early
stages of the scientific process.

### Contributions

We present ‘cyber-physicist’ (CyPhy), our novel AI research associate for
early-stage scientific process of hypothesis generation and initial
(in)validation, grounded in the most invariable mathematical foundations of
classical and relativistic physics. Our framework distinguishes itself from
existing rule-based reasoning, statistical learning, and hybrid AI methods by:

  1. (1)

an ability to rapidly enumerate and test a diverse set of mathematically sound
and parsimonious physical hypotheses, starting from a few basic assumptions on
the embedding spacetime topology;

  2. (2)

a distinction between non-negotiable mathematical truism (e.g., conservation
laws or symmetries), that are directly implied by properties of spacetime, and
phenomenological relations (e.g., constitutive laws), whose characterization
relies indisputably on empirical observation, justifying targeted use of data-
driven methods (e.g., ML or polynomial regression); and

  3. (3)

a “simple-first” strategy (following Occam’s razor) to search for new
hypotheses by incrementally introducing latent variables that are expected to
exist based on topological foundations of physics.

### Background

AI-assisted discovery of scientific knowledge has been an active area of
research (Langley 1998) long before the rise of GPU-accelerated deep learning
(DL). As computational power and data sources are becoming more ubiquitous,
model-based, data-driven, and hybrid AI methods are playing an increasingly
more important role in various scientific activities (Kitano 2016; Raghu and
Schmidt 2020).

Related efforts to our approach to scientific hypothesis generation and
evaluation are mostly engineered after how humans approach scientific
discovery, including sequential rule-based symbolic regression (Schmidt and
Lipson 2009; Udrescu and Tegmark 2020), latent space representation learning
via deep neural net auto-encoders (Iten et al. 2020; Nautrup et al. 2020) and
strategic combinations of divide-and-conquer, unsupervised learning,
simplification by penalizing description lengths in the loss function, and a
posteriori unification by clustering (Wu and Tegmark 2019). While these and
other efforts have shown great promise for elevating AI to the role of an
autonomous, creative, and insightful collaborator that can offer human
scientists a set of viable options to consider, their applications have
remained limited to rather basic examples.

On the more domain-specialized end, DL has been widely successful in
classification, regression, and forecasting tasks in scientific areas as
diverse as turbulence (Miyanawala and Jaiman 2017; Wang et al. 2020), chaotic
particle dynamics (Breen et al. 2020), molecular chemistry and materials
science (Butler et al. 2018), and protein engineering (Yang, Wu, and Arnold
2019), among others. Most specialized DL architectures are ad hoc, designed
(by humans) using narrow, domain-specific, and (by construction) biased
knowledge and expertise, stifling innovation and surprise. Moreover, DL models
that successfully capture nontrivial patterns in data are often difficult to
explain, lack guarantees even within their training space, and poorly
extrapolate to out-of-training scenarios (Mehta et al. 2019). Training such
models for high-dimensional physics problems requires enormous data, which is
either unavailable or too costly to obtain in many experimental sciences.

## The Cyber-Physicist

We introduce an AI tool that can bridge multiple levels of abstraction, using
a domain-agnostic representation scheme to express a wide range of
mathematically viable physical hypotheses by exploiting common structural
invariants across physics. Our approach entails:

  1. (a)

defining a relatively unbiased ontology rooted in fundamental abstractions
that are common to all known theories of classical and relativistic physics;

  2. (b)

constructing a constrained search space to enumerate viable hypotheses with
postulated invariants, e.g., built-in conservation laws that are consistent
with the presupposed spacetime topology; and

  3. (c)

automatically assembling interpretable ML architectures for each hypothesis,
to estimate parameters for phenomenological relations from empirical data.

At the core of (a) is a powerful mathematical abstraction of physical
governing equations rooted in algebraic topology and differential geometry
(Frankel 2011). This abstraction leads to an ontological commitment to the
relationship between physical measurement and basic properties of the
embedding spacetime—but nothing more, to leave room for innovation and
surprise. This relationship has been shown to be responsible for the analogies
and common structure across physics (Tonti 2013), exploited in (b), along with
search heuristics based on analogical reasoning. Each viable hypothesis is
automatically compiled to an interpretable “computation graph”—tensor-based
architecture, akin to a neural net with convolution layers to compute
differentiation/integration and (non)linear local operators for constitutive
equations—for a given cellular decomposition of embedding spacetime using
well-established concepts from cellular homology (Hatcher 2001) and exterior
calculus of differential and discrete forms (Bott and Tu 1982; Hirani 2003)
that are under-utilized in AI.

### Topological Foundations of Physics

The key enabler of our AI framework is a simple type system for (a) physical
variables, based on how they are measured in spacetime; and (b) physical
relations, based on their (topological vs. metric) nature, and the variables
they connect. Following the ground-breaking discoveries by a number of
mathematicians, physicists, and electrical engineers (Kron 1963; Roth 1955;
Branin 1966) towards a general network theory, Tonti explained the fascinating
analogies across classical and relativistic physics in his pioneering life-
long work (Tonti 2013) by reframing them in the language of cellular homology,
leading to informal classification diagrams. Tonti diagrams can be formalized
as directed graphs with strongly typed nodes for variables and edges for
relations. The variable are typed as $(\mathrm{d}_{1},\mathrm{d}_{2})-$forms
based on their measurement on $\mathrm{d}_{1}-$ and
$\mathrm{d}_{2}-$dimensional submanifolds ($\mathrm{d}_{1}-$ and
$\mathrm{d}_{2}-$cells) of space and time, respectively. For instance, to
model heat transfer in (3+1)D spacetime, temperature is typed as a
$(0,1)-$form because it is measured at spatial points ($0-$cells) and during
temporal intervals ($1-$cells), whereas heat flux is typed as a $(2,1)-$form
because it is measured over spatial surfaces ($2-$cells) and during temporal
intervals ($1-$cells). In classical calculus, both of these variables reduce
to scalar and vector fields, probed at spatial points and at temporal
instants, to write down compact pointwise DEs; however, keeping track of the
topological and geometric character of DEs is key to a deeper understanding of
how known physical theories work, and building on top of it for AI-assisted
discovery of new physics grounded in mathematical foundations.

Figure 1: A topology-aware representation for physics (Tonti 2013): (a)
variables associated with spatial and temporal cells of various dimensions
give rise to primary forms and secondary forms (also called pseudo-forms); (b)
resulting in 32 possible types for spatio-temporal forms, and an underlying
structure for fundamental theories of physics.

The spatiotemporal cells (or embedding manifolds) are further classified as
primary or secondary, endowed with inner or outer orientations, respectively,
depending on how the variables change sign in a hypothetical reversal of
spacetime orientation (Mattiussi 2000). The cells are related by topological
duality (Fig. 1 (a)). For example, an inner-oriented curve ($1-$cell,
$\overline{\sigma}^{1}$) sitting in primary space, along which temperature
variations are measured, is dual to an outer-oriented surface ($2-$cell,
$\widetilde{\sigma}^{2}$) sitting in secondary space, over which heat flux is
measured, and the two cells are spatially registered and consistently
oriented, if we embed them in two co-located “copies” of 3D space.

Figure 2: Tonti diagrams are recipes to generate governing equations in
different contexts, defined by a continuum, discrete, or semi-discrete setting
and a topological embedding of the variables based on how they are measured.
The conservation laws in terms of co-boundary operators result directly from
assumed properties of space (or spacetime), while constitutive relations must
be learned from data (e.g., via regression/ML).

The relations among variables on the Tonti diagrams are typed based on the
pairs of variables they relate, as well as the nature of the relation itself:

  * •

Topological relations map spatiotemporal forms to forms of one higher
dimension in space or time via incidence relations, and are responsible for
propagation of information in spacetime through incident cells.

  * •

Metric relations locally map forms defined over dual cells to one another
based on phenomenological properties, spatial lengths, and temporal durations,
and are responsible for local distortion of information.

  * •

Algebraic relations are in-place, i.e., map a given form to another from of
the same type, and can be used to capture initial/boundary conditions,
external source/sink terms, or cross-physics couplings between variables of
the same type on different diagrams.

The relations are drawn in Fig. 1 as vertical arrows, horizontal (or
horizontal-diagonal) arrows, and loops ($1-$cycles), respectively. The
interpretation of these relations to symbolic or numerical operations depends
on the choice of a cellular decomposition of spacetime on which they operate.
For example, using a continuum spacetime with infinitesimal cells, the
variables are viewed as differential forms and the topological operators on
them are interpreted as exterior derivatives (Bott and Tu 1982). In elementary
calculus, these operators give rise to gradient, curl, and divergence in space
and partial derivative in time in terms of scalar and vector fields that are
proxy to these forms, leading to partial DEs (PDEs). In a discrete (or semi-
discrete) setting, on the other hand, the same diagram can be used to produce
integral (or integro-differential) equations that capture the same fundamental
conservation and constitutive realities, where the variables are viewed as co-
chains, also called discrete forms (or mixed forms, e.g., discrete in space,
differential in time, or vice versa) and the topological operators become co-
boundary operators that are fundamental in cellular homology (Hatcher 2001).
For example, using a semi-discretization in space with integral quantities
associated with $0-$, $1-$, $2-$ and $3-$cells on a pair of staggered
unstructured meshes in 3D, while keeping time as a continuum, the semi-
discrete form of the heat equation as a system of ordinary DEs (ODEs) (Fig. 2
(a)). Upon discretization of time, one obtains algebraic equations that can be
solved or parameter estimated via tensor-based ML.

Figure 3: Tonti diagrams capture the common structure responsible for
analogies across classical and relativistic physics with a clear distinction
between topological and phenomenological relations that follow certain rules.

It is important to note that 3D meshes in space and 1D time-stepping are not
the only ways to provide a combinatorial topology to interpret Tonti diagrams
in a discrete setting. Another example is a directed graph representation of
lumped-parameter networks such as system models in Modelica or electrical
circuits in Spice. The variables in this case are associated with nodes,
edges, and meshes (i.e., primitive cycles) and incidence relations are
obtained from graph connectivity and edge directions. The same topological
operator that leads to a spatial divergence, discretized by a sum of fluxes on
the incident faces of a volume in a 3D mesh, also leads to a superposition of
forces on interacting planets, sum of currents in/out of junctions in
electrical circuits, and superposition of torques on kinematic chains (Fig. 2
(b, c, d)). Both ODEs and PDEs and their integral or integro-differential
forms upon full or semi-discretization can be captured with the same
(abstract) operators, and Tonti diagrams serve as recipes to compose them to
generate governing equations.

Figure 3 shows a few other examples of Tonti diagrams for fundamental theories
in classical and relativistic physics. The differences amount to (a)
topological and metric context; (b) relevant variables and their
dimensions/units; and (c) phenomenological relations.

Figure 4: I-nets are generalizations of Tonti diagrams for finite topological
products of finite-dimensional spaces with relaxed rules for feasible
phenomenological links to accomodate middle-ground theories.

### An Ontology for Scientific Process

We present a novel representation, called ‘interaction networks’ (I-nets),
based on a generalization of Tonti diagrams that is expressive and versatile
enough to accommodate novel scientific hypotheses, while retaining a basic
commitment to philosophical principles such as parsimony (Occam’s razor),
measurement-driven classification of variables, and separation of non-
negotiable mathematical properties of spacetime (homology) from domain-
specific empirical knowledge (phenomenology). Data science is employed to help
only with the latter.

  * •

We conceptualize three levels of abstraction related by inheritance: abstract
(symbolic) I-nets $\to$ discrete (cellular) I-nets $\to$ numerical (tensor-
based) I-nets.

  * •

At each level, an I-net instance is contextualized by user-defined assumptions
on spacetime topology, semantics of physical quantities, and structural
restrictions on allowable diagrams based on analogical reasoning and domain-
specific insight (if available).

  * •

Every I-net instance distinguishes between topological and metric operators;
however, it has additional degrees of freedom (beyond Tonti diagrams) for the
latter to allow for phenomenological relations among variables that may not be
dual to each other.

The latter is motivated by the observation that some existing middle-ground
theories use phenomenological relations to capture a combination of
topological and metric aspects.

We define an abstract (symbolic) I-net on a single $\mathscr{D}-$space as a
finite collection of primary and/or secondary co-chain complexes that are
inter-connected by phenomenological links, as shown in Fig. 4 (a). Each co-
chain complex is a sequence of (symbolic) $\mathrm{d}-$forms related by
(symbolic) co-boundary operators from $\mathrm{d}-$forms to
$(\mathrm{d}+1)-$forms ($0\leq\mathrm{d}\leq\mathscr{D}$). The interpretation
of $\mathrm{d}\to(\mathrm{d}+1)$ maps depends on the embedding dimension
$\mathscr{D}$; for instance, if $\mathscr{D}=1$ the only option for the input
is $\mathrm{d}=0$ leading to a simple partial derivative ($0\to 1$), whereas
for $\mathscr{D}=3$, we can have $\mathrm{d}=0,1,2$ leading to gradient ($0\to
1$), curl ($1\to 2$), and divergence ($2\to 3$) operations, respectively.

These sequences may represent different (mechanical, electrical, thermal,
etc.) domains of physics. Although, for most known physics, each domain’s
theory appears as one pair of (primary and secondary) sequences in tandem,
connected by horizontal (or horizontal-diagonal) constitutive relations
leading to Tonti diagrams, we do not make any such restriction when looking
for new theories. The cross-sequence links can thus represent both single-
physics constitutive relations and mutli-physics coupling interactions.
Conservation laws, on the other hand, are represented by a balance between the
output of a topological operator and an external source/sink, the latter being
represented by a loop.

It is often more convenient to define product spaces (e.g., separate 3D space
and 1D time, as opposed to 4D spacetime) in which conservation laws are stated
as sums of incoming topological relations being balanced against an external
source/sink. To accommodate such representations, we define abstract
(symbolic) I-nets on a product of a $\mathscr{D}_{1}-$space and a
$\mathscr{D}_{2}-$space as multi-sequences of co-chains, connected by
phenomenological links, as before. It is possible to form $2^{2}=4$ possible
such multi-sequences with various orientation combinations, two of which lead
to so-called mechanical and field theories (Tonti 2013), shown in Fig. 1 for
(3+1)D spacetime and repeated in Fig. 4 (b) for higher-dimensional pairs of
abstract topological spaces. This construction is generalized to products of
more than two spaces in a straightforward combinatorial fashion.

Based on the topological context, the semantics for co-boundary operators is
unambiguously determined by the dimensions of the two variables (i.e., co-
chains) they relate. However, phenomenological links require specifying a
parameterization of possibly nonlinear, in-place, and purely metric relations
they represent, using unknown parameters that must be learned from data.

Figure 5: The search space for the dynamics of a pendulum in 1D time. The
complete hypotheses (yellow nodes) correspond to I-net structures that pose
new nontrivial equations to be tested against data, whereas incomplete
hypotheses (white nodes) have “dangling” branches that are completed in their
child states.

Once one or more hypotheses are specified in the language of abstract
(symbolic) I-nets with unknown phenomenological parameters (e.g., thermal
conductivity in the earlier heat transfer example), the parameters can be
optimized to fit the data and the regression error can be used to evaluate the
fitness of hypotheses.

### A Search for Viable Hypotheses

Having defined a combinatorial representation of viable hypotheses that are
partially ordered in terms of complexity, the next step is to generate and
test the hypotheses in a “simple-first” fashion. The search space is defined
by a directed acyclic graph (DAG) whose nodes (i.e., ‘states’) represent
symbolic I-net instances. The edges (i.e., state transitions) represent
generating a new I-net structure by incrementally adding complexity to the
parent state. Each action can be one or composition of (a) defining a new
symbolic variable, in an existing co-chain complex, by applying a topological
operator to an existing variable; (b) defining a new variable in a latent co-
chain complex; and (c) adding phenomenological links of prescribed form and
unknown parameters, connecting existing variables. The search is guided by a
loss function determined by how well the hypotheses represented by these I-net
structures explain a given dataset. The algorithm may also be equipped with
user-specified heuristic rules to prune the search space or prioritize paths
that are perceived as “more likely” due to structural analogies with existing
theories.

Figure 6: The hypotheses H-04 and H-08 of Fig. 5 are enumerated and visualized
by the software and evaluated against data (split 0.7-0.3 for
training/testing). Both energy (first-order) and and torque (second-order)
forms of the governing equation are discovered without human intervention. The
former was quite unexpected, since its I-net structure does not correspond to
a Tonti diagram. The latter has a larger error due to finite difference
discretization.

The input to the search algorithm includes the bare minimum contextual
information such as the assumed underlying topology, a preset number of
physical domains, and the types of measured variables, e.g., spatiotemporal
associations, tensor ranks and shapes, and dimensions/units. The search starts
from an “initial” I-net instance (i.e., the ‘root’) that embodies only
measured variable(s) with no initial edges except the ones that are asserted a
priori, e.g., loops for initial/boundary conditions or source terms, if
applicable. The spatio-temporal types and physical semantics for these
variables are provided by the experimentalist.

For example, consider a simple pendulum (Fig. 5 (a)). We have only 1D time,
leading to a topological space of inter-connected time instants
$\overline{\tau}^{0},\widetilde{\tau}^{0}=\overline{\tau}^{0}+\nicefrac{{\epsilon}}{{2}}$
and time intervals
$\overline{\tau}^{1}=(\overline{\tau}^{0},\overline{\tau}^{0}+\epsilon),\widetilde{\tau}^{1}=(\widetilde{\tau}^{1},\widetilde{\tau}^{1}+\epsilon)$
to which data may be associated. Suppose we are given time series data for
angular position $\theta(\overline{\tau}^{0})$. The initial I-net instance is
a single symbolic variable for this $0-$form, which can be differentiated only
once in primary 1D time to obtain angular velocity as a $1-$form:
$\theta(\overline{\tau}^{0})\to\omega(\overline{\tau}^{1})=\delta[\theta](\overline{\tau}^{1})$
at the root of the search DAG (Fig. 5 (b)). The DAG is expanded by adding new
phenomenological links, either between two existing variables, or between an
existing variable and one in a newly added latent co-chain sequence (Fig. 5
(c)). In this example, the hypotheses are numbered H-00 (the root) through
H-15, enumerating all possible I-net structures formed by at most one latent
co-chain complex in 1D time. The user can specify the maximum number of latent
variables that the algorithm may consider, to keep the search tractable.

Not every introduction of new variables or relations makes nontrivial
statements about physics. For example, the hypothesis H-01 produces a new
variable typed as a $1-$pseudo-form
$T(\widetilde{\tau}^{1})=\mathcal{f}_{1}(\theta(\ast\widetilde{\tau}^{1}))$,
where the $\ast-$operator takes $\widetilde{\tau}^{1}$ to its dual:
$\ast(\widetilde{\tau}^{0},\widetilde{\tau}^{0}+\epsilon)=\widetilde{\tau}^{0}+\nicefrac{{\epsilon}}{{2}}$.
However, until this new variable is reached through another path to close a
cycle and pose a nontrivial equation, we do not have a complete hypothesis to
(in)validate against data. Further down the search DAG, H-08 defines a new
variable typed as a $0-$pseudo-form
$L(\widetilde{\tau}^{0})=\mathcal{f}_{2}(\omega(\ast\widetilde{\tau}^{0}))$
where
$\ast\widetilde{\tau}^{0}=(\widetilde{\tau}^{0}-\nicefrac{{\epsilon}}{{2}},\widetilde{\tau}^{0}+\nicefrac{{\epsilon}}{{2}})$.
The co-boundary operation $L(\widetilde{\tau}^{0})\to
T(\widetilde{\tau}^{1})=\delta[L](\widetilde{\tau}^{1})$, closes the cycle and
produces a commutative diagram (Fig. 5 (c)) leading to:

$\mathcal{E}_{\textsf{H-08}}(\theta;\mathcal{f}_{1},\mathcal{f}_{2})=\mathcal{f}_{1}(\theta)-\delta^{\ast}[\mathcal{f}_{2}(\delta[\theta])]=0,$
(1)

where $\mathcal{f}_{1},\mathcal{f}_{2}$ are selected from restricted function
spaces $\mathcal{F}_{1},\mathcal{F}_{2}$ to avoid overfitting (e.g.,
parameterized by a linear combination of domain-aware basis functions) and
their parameters must be determined from data to minimize the residual error
$\mathcal{E}_{\textsf{H-08}}$ over the entire period of data collection. A
loss function can, for example, be defined as a mean-squared-error (MSE) to
penalize violations uniformly over the time series period:

$\textsf{Loss}_{\textsf{H-08}}=\min_{\mathcal{f}_{1}\in\mathcal{F}_{1}}\min_{\mathcal{f}_{2}\in\mathcal{F}_{2}}\big{\|}\mathcal{E}_{\textsf{H-08}}(\theta;\mathcal{f}_{1},\mathcal{f}_{2})\big{\|}_{\widetilde{\tau}^{1}},$
(2)

where $\|\cdot\|_{\widetilde{\tau}^{1}}$ is an $L_{2}-$norm computed as a
temporal integral, i.e., sum of squared errors
$\mathcal{E}_{\textsf{H-08}}^{2}(\theta;\mathcal{f}_{1},\mathcal{f}_{2})$ over
time intervals $\widetilde{\tau}^{1}$ where (1) is evaluated. In this example,
it turns out that the best fit is achieved with
$\mathcal{f}_{1}(\theta)=\mathcal{c}_{1}\sin\theta$ and
$\mathcal{f}_{2}(\omega)=\mathcal{c}_{2}\omega$ where
$\nicefrac{{\mathcal{c}_{2}}}{{\mathcal{c}_{1}}}=-\nicefrac{{g}}{{r}}$. The
latent variables $L(\widetilde{\tau}^{0})$ and $L(\widetilde{\tau}^{0})$ turn
out to be the familiar notions of angular momentum and torque, respectively,
although the software need not know anything about them to generate and test
what-if scenarios about their existence and correlations with angular position
and velocity. Hence, interpretability of the discovered relationships by a
human scientist does not require predisposing the AI associate to such
interpretations, enabling unexpected discoveries.

In general, every state in the search DAG can be classified as complete or
incomplete hypotheses. The former are I-net structures with “dangling”
branches that carry no new nontrivial information in addition to their parent
states. Every time such a branch is turned into one or more closed cycles by
adding enough new variables and/or relations, a new constraint is hypothesized
that can be evaluated against data. When adding new dangling branches to the
I-net structure, the search algorithm prioritizes actions that produce I-net
structures similar to existing Tonti diagrams by assigning a penalty factor to
every violation of the common structure (e.g., diagonal phenomenological links
connecting non-dual cells). The loss for complete hypotheses can be computed
as the sum of penalties for the I-net structure and the sum of residual errors
for each of the independent constraints, implied by converging paths,
multiplied by use-specified relative weight of the penalties and errors. We
use an A* algorithm to search the space of hypotheses. Since we cannot compute
the error for incomplete hypotheses, we can only prune them when the increase
in their penalty is large enough that it would fail even if it had no error at
all.

### Generating Symbolic Expressions

One of the practical features of our implementation in Python is its ability
to automatically convert I-net instances to symbolic DE expressions in SymPy,
when the co-boundary operators are interpreted in a differential setting for
infinitesimal cells ($\epsilon\to 0^{+}$); for example, equation (1) can be
rewritten as a nonlinear ODE:

$\mathcal{E}_{\textsf{H-08}}(\theta;\mathcal{f}_{1},\mathcal{f}_{2})=\mathcal{f}_{1}(\theta)-\frac{\partial}{\partial
t}\left[\mathcal{f}_{2}\left(\dot{\theta}(t)\right)\right].$ (3)

As a result, the generated hypotheses can be evaluated using any number of
existing ML or symbolic regression frameworks that standardize on ODE/PDE
inputs. For example, using non-orthogonal basis functions $\\{1,x,x^{2},\sin
x,\cos x\\}$ to span both function spaces $\mathcal{F}_{1},\mathcal{F}_{2}$,
we can substitute for both symbolic functions:

$\displaystyle\mathcal{f}_{1}(\theta)$
$\displaystyle:=\mathcal{c}_{0}^{1}+\mathcal{c}_{1}^{1}\theta+\mathcal{c}_{2}^{1}\theta^{2}+\mathcal{c}_{3}^{1}\sin\theta+\mathcal{c}_{3}^{1}\cos\theta,$
(4) $\displaystyle\mathcal{f}_{2}(\dot{\theta})$
$\displaystyle:=\mathcal{c}_{0}^{2}+\mathcal{c}_{1}^{2}\dot{\theta}+\mathcal{c}_{2}^{2}\dot{\theta}^{2}+\mathcal{c}_{3}^{2}\sin\dot{\theta}+\mathcal{c}_{3}^{2}\cos\dot{\theta},$
(5)

into (3) to obtain a symbolic second-order (non)linear ODE in SymPy. Next, the
software performs algebraic simplification to identify equivalence classes of
hypotheses that, despite coming from different I-net structures, lead to the
same ODE upon differential interpretation of the I-nets. For ODEs which, after
simplification, are linear combinations of nonlinear (differential/algebraic)
terms that are computable from data, we can apply symbolic regression to
estimate the coefficients from data; for example, we tried LASSO-regularized
least-squares regression in PDE-FIND (Rudy et al. 2017) where each term
involving a derivative is evaluated using finite difference or polynomial
approximation, whose results are shown in Fig. 6.

There are at least two issues with this approach:

First, more sophisticated regression or nonlinear programming methods are
needed if the DE has terms that have nested nonlinear functions, i.e., cannot
be represented as a linear combination of nonlinear terms because of unknown
coefficients embedded within each term. We solve this problem by directly
mapping I-net structures to computation graphs in PyTorch, skipping
differential interpretation to symbolic DEs altogether.

Second, numerical approximation of symbolic PDEs is a tricky business, as the
discrete forms (in 3D space) may not obey the conservation principles
postulated by the I-net structure after such approximations. It is difficult
to separate discretization errors from modeling errors and noise in data. One
of the key advantages of I-nets is the rich geometric information in their
type system that is fundamental to physics-compatible and mimetic
discretization schemes (Koren et al. 2014; Palha et al. 2014; Lipnikov,
Manzini, and Shashkov 2014) that ensure conservation laws are satisfied
exactly as a discrete level, regardless of spatial mesh or time-step
resolutions. Such information is lost upon conversion to symbolic DEs.
Retaining this information is even more important when dealing with noisy
data, because discrete differentiation of noisy data (e.g., via finite
difference or polynomial fitting) can substantially amplify the noise.

The good news is that we can directly interpret the same I-net instance in
integral form to generate equations over larger regions in space and/or time,
to make the computations more resilient to noise. For example, in the heat
equation, the discrete divergence of heat flux over a single $3-$cell is
replaced by a flux integral over a collection of $3-$cells, and is equated
against the volumetric intgeral of internal energy within the collection. The
cancellation of internal surface fluxes (discrete form of Gauss’ divergence
theorem) is built into the interpretation based on cellular homology. The
integrals can be computed using higher-order integration schemes, e.g., using
polynomial interpolation with underfitting to filter the noise.

Further details on directly and automatically mapping the abstract (symbolic)
I-net structures to discrete (cellular) and numerical (tensor-based) I-net
instance (e.g., computation graphs in PyTorch), learning scale-aware
phenomenological relations, and physics-compatible discretization and de-
noising will be presented in a full paper.

Figure 7: The search DAG and a number of viable hypotheses to explain
ultrasound wavefield in metal parts. Figure 8: The AI associate discovers the
(integral form) of the wave equation as well as the proper length/time scale
at which the heterogeneous material properties (in this case, speed of sound)
must be defined.

## Real-World Scientific Discovery

Figures 7 and 8 illustrate the application of our AI approach to an
elastodynamics challenge problem provided by AFRL in the course of the DARPA
AI Research Associate (AIRA) program that supported the development of CyPhy.
The input is noisy data obtained by ultrasound imaging, measured in (2+1)D
spacetime over the surface of several material samples with heterogeneous
properties.

Figure 7 illustrates the search DAG along with a number of I-net structures
for viable hypotheses, each postulating the relevance of a conservation law
and existence of a few phenomenological relations. Figure 8 shows the ranking
of these hypotheses based on their residual errors when tested against data.
Each hypothesis can be interpreted in differential, integral, or integro-
differential forms. The results demonstrate that integral forms applied to
wide spatial and temporal neighborhoods (of $\sim 25$ grid elements along each
axis) with high-order polynomial underfitting (up to cubic in each
coordinate), resulting in a length/time scale-aware definition of (nonlocal)
phenomenological relations as well as physics-compatible (i.e., mimetic)
discretization and de-noising, are preferable to strictly local numerical
schemes such as finite difference discretization.

## Conclusion

Statistical learning methods, despite their accuracy and efficiency in narrow
regimes for which they are carefully engineered, are not sufficient to
independently acquire deep understandings of the scientific problems they are
applied to. Human scientists continue to handle most of knowledge-centric
aspects of the scientific process based on domain-specific insight,
experience, and expertise.

Our novel approach to early-stage scientific hypothesis generation and testing
demonstrates a path forward towards context-aware, generalizable, and
interpretable AI for scientific discovery. Our AI associate (CyPhy)
distinguishes between non-negotiable mathematical truism, implied by the
relationship between measurement and presupposed spacetime topology, and
phenomenological realities that are at the mercy of empirical learning. Data-
driven regression is targeted at the latter to enable distilling governing
equations from sparse and noisy data, while providing deep insights into the
mathematical foundations.

## Acknowledgment

This material is based upon work supported by the Defense Advanced Research
Projects Agency (DARPA) under Agreement No. HR00111990029.

## References

  * Bott and Tu (1982) Bott, R.; and Tu, L. W. 1982.  _Differential Forms in Algebraic Topology_.  Springer Science & Business Media.
  * Branin (1966) Branin, F. H. 1966.  The Algebraic-Topological Basis for Network Analogies and the Vector Calculus.  In _Symposium on Generalized Networks_ , 453–491. Polytechnic Institute of Brooklyn, NY.
  * Breen et al. (2020) Breen, P. G.; Foley, C. N.; Boekholt, T.; and Zwart, S. P. 2020.  Newton versus the Machine: Solving the Chaotic Three-body Problem Using Deep Neural Networks.  _Monthly Notices of the Royal Astronomical Society_ 494(2): 2465–2470.
  * Butler et al. (2018) Butler, K. T.; Davies, D. W.; Cartwright, H.; Isayev, O.; and Walsh, A. 2018.  Machine Learning for Molecular and Materials Science.  _Nature_ 559(7715): 547–555.  doi:10.1038/s41586-018-0337-2.
  * Cranmer et al. (2019) Cranmer, M. D.; Xu, R.; Battaglia, P.; and Ho, S. 2019.  Learning Symbolic Physics with Graph Networks.  _arXiv preprint arXiv:1909.05862_ .
  * Daw et al. (2020) Daw, A.; Thomas, R. Q.; Carey, C. C.; Read, J. S.; Appling, A. P.; and Karpatne, A. 2020.  Physics-Guided Architecture (PGA) of Neural Networks for Quantifying Uncertainty in Lake Temperature Modeling.  In _Proceedings of the 2020 SIAM International Conference on Data Mining_ , 532–540. Society for Industrial and Applied Mathematics (SIAM).
  * Frankel (2011) Frankel, T. 2011.  _The Geometry of Physics: An Introduction_.  Cambridge University Press.
  * Hatcher (2001) Hatcher, A. 2001.  _Algebraic Topology_.  Cornell University.
  * Hirani (2003) Hirani, A. N. 2003.  _Discrete Exterior Calculus_.  Ph.D. dissertation, California Institute of Technology.
  * Iten et al. (2020) Iten, R.; Metger, T.; Wilming, H.; Del Rio, L.; and Renner, R. 2020.  Discovering Physical Concepts with Neural Networks.  _Physical Review Letters_ 124(1): 010508.
  * Kitano (2016) Kitano, H. 2016.  Artificial Intelligence to Win the Nobel Prize and Beyond: Creating the Engine for Scientific Discovery.  _AI Magazine_ 37(1): 39–49.
  * Koren et al. (2014) Koren, B.; Abgrall, R.; Bochev, P.; Frank, J. E.; and Perot, B. 2014.  Physics-Compatible Numerical Methods.  _Journal of Computational Physics_ 257(Part B): 1039–1039.
  * Kron (1963) Kron, G. 1963.  _Diakoptics: The Piecewise Solution of Large-Scale Systems_ , volume 2.  MacDonald.
  * Langley (1998) Langley, P. 1998.  The Computer-Aided Discovery of Scientific Knowledge.  In _International Conference on Discovery Science_ , 25–39. Springer.
  * Lipnikov, Manzini, and Shashkov (2014) Lipnikov, K.; Manzini, G.; and Shashkov, M. 2014.  Mimetic Finite Difference Method.  _Journal of Computational Physics_ 257: 1163–1227.
  * Mattiussi (2000) Mattiussi, C. 2000.  The Finite Volume, Finite Element, and Finite Difference Methods as Numerical Methods for Physical Field Problems.  _Advances in Imaging and Electron Physics_ 113: 1–147.
  * Mehta et al. (2019) Mehta, P.; Bukov, M.; Wang, C. H.; Day, A. G. R.; Richardson, C.; Fisher, C. K.; and Schwab, D. J. 2019.  A High-Bias, Low-Variance Introduction to Machine Learning for Physicists.  _Physics Reports_ 810: 1–124.
  * Miyanawala and Jaiman (2017) Miyanawala, T. P.; and Jaiman, R. K. 2017.  An Efficient Deep Learning Technique for the Navier-Stokes Equations: Application to Unsteady Wake Flow Dynamics.  _arXiv Preprints_ ISSN 0264-6021.  doi:10.1016/j.eswa.2008.08.077.
  * Nautrup et al. (2020) Nautrup, H. P.; Metger, T.; Iten, R.; Jerbi, S.; Trenkwalder, L. M.; Wilming, H.; Briegel, H. J.; and Renner, R. 2020.  Operationally Meaningful Representations of Physical Systems in Neural Networks.  _arXiv preprint arXiv:2001.00593_ .
  * Palha et al. (2014) Palha, A.; Rebelo, P. P.; Hiemstra, R.; Kreeft, J.; and Gerritsma, M. 2014.  Physics-Compatible Discretization Techniques on Single and Dual Grids, with Application to the Poisson Equation of Volume Forms.  _Journal of Computational Physics_ 257: 1394–1422.
  * Raghu and Schmidt (2020) Raghu, M.; and Schmidt, E. 2020.  A Survey of Deep Learning for Scientific Discovery.  _arXiv preprint arXiv:2003.11755_ .
  * Raissi, Perdikaris, and Karniadakis (2019) Raissi, M.; Perdikaris, P.; and Karniadakis, G. E. 2019.  Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations.  _Journal of Computational Physics_ 378: 686–707.
  * Roth (1955) Roth, J. P. 1955.  An Application of Algebraic Topology to Numerical Analysis: On the Existence of a Solution to the Network Problem.  _Proceedings of the National Academy of Sciences_ 41(7): 518–521.
  * Rudy et al. (2017) Rudy, S. H.; Brunton, S. L.; Proctor, J. L.; and Kutz, J. N. 2017.  Data-Driven Discovery of Partial Differential Equations.  _Science Advances_ 3(4): e1602614.
  * Sanchez-Gonzalez et al. (2020) Sanchez-Gonzalez, A.; Godwin, J.; Pfaff, T.; Ying, R.; Leskovec, J.; and Battaglia, P. W. 2020.  Learning to Simulate Complex Physics with Graph Networks.  _arXiv preprint arXiv:2002.09405_ .
  * Schmidt and Lipson (2009) Schmidt, M.; and Lipson, H. 2009.  Distilling Free-Form Natural Laws from Experimental Data.  _Science_ 324(5923): 81–85.
  * Seo and Liu (2019) Seo, S.; and Liu, Y. 2019.  Differentiable Physics-Informed Graph Networks.  _arXiv preprint arXiv:1902.02950_ .
  * Stevens et al. (2020) Stevens, R.; Taylor, V.; Nichols, J.; Maccabe, A. B.; Yelick, K.; and Brown, D. 2020\.  AI for Science.  Technical report, Argonne National Laboratory(ANL).
  * Tonti (2013) Tonti, E. 2013.  _The Mathematical Structure of Classical and Relativistic Physics: A General Classification Diagram_.  Modeling and Simulation in Science, Engineering and Technology. Birkhäuser.  ISBN 9781461474210.
  * Udrescu and Tegmark (2020) Udrescu, S. M.; and Tegmark, M. 2020.  AI Feynman: A Physics-Inspired Method for Symbolic Regression.  _Science Advances_ 6(16): eaay2631.
  * Wang et al. (2020) Wang, R.; Kashinath, K.; Mustafa, M.; Albert, A.; and Yu, R. 2020.  Towards Physics-Informed Deep Learning for Turbulent Flow Prediction.  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 1457–1466.
  * Wei and Chen (2019) Wei, Z.; and Chen, X. 2019.  Physics-Inspired Convolutional Neural Network for Solving Full-Wave Inverse Scattering Problems.  _IEEE Transactions on Antennas and Propagation_ 67(9): 6138–6148.
  * Wu and Tegmark (2019) Wu, T.; and Tegmark, M. 2019.  Toward an AI Physicist for Unsupervised Learning.  _arXiv preprint arXiv:1810.10525_ .
  * Yang, Wu, and Arnold (2019) Yang, K. K.; Wu, Z.; and Arnold, F. H. 2019.  Machine-Learning-Guided Directed Evolution for Protein Engineering.  _Nature Methods_ 16(8): 687–694.  doi:10.1038/s41592-019-0496-6.


#  Unconventional rheological properties in systems of deformable particles

Anshuman Pasupalak  Shawn Khuhan Samidurai  Yanwei Li  Yuanjian Zheng
Division of Physics and Applied Physics, School of Physical and Mathematical
Sciences
Nanyang Technological University  Ran Ni  Division of Physics and Applied
Physics, School of Physical and Mathematical Sciences
Nanyang Technological University  School of Chemical and Biomedical
Engineering
Nanyang Technological University  Massimo Pica Ciamarra  Division of Physics
and Applied Physics, School of Physical and Mathematical Sciences
Nanyang Technological University  CNR–SPIN, Dipartimento di Scienze Fisiche,
Università di Napoli Federico II, I-80126, Napoli, Italy

###### Abstract

We demonstrate the existence of unconventional rheological and memory
properties in systems of soft-deformable particles whose energy depends on
their shape, via numerical simulations. At large strains, these systems
experience an unconventional shear weakening transition characterized by an
increase in the mechanical energy and a drastic drop in shear stress, which
stems from the emergence of short-ranged tetratic order. In these weakened
states, the contact network evolves reversibly under strain reversal, keeping
memory of its initial state, while the microscopic dynamics is irreversible.

##  I Introduction

When compressed to high density, soft particles such as microgels, foams and
cells in tissue deform, develop facets and tessellate available space,
provided they do not interpenetrate. In these high-density systems, pair-wise
interacting descriptions lose their value [1] as the energy of a particle
depends on the positions of all its neighbors and is defined through a many-
body potential. This potential is the work needed to squeeze the particle into
its volume, and establishes a close link between the mechanical and the
geometrical properties of the system. This link has recently been uncovered in
a two-dimensional model, when the energy of a particle depends on its
perimeter [2, 3, 4, 5, 6]. In this case, an elementary relaxation event
involving the disappearance of an edge separating two particles of length $l$,
has an energy cost $\propto l$. Hence, the system’s mechanical response
correlates with the edge-length distribution, $P(l)$, and in particular, the
system loses mechanical rigidity when an extensive number of zero-length edges
appear, $P(l)\propto l^{\theta}$ with $\theta=0$ [2, 3, 4, 5, 6]. When a
particle’s energy depends not only on its perimeter but also on other
geometrical properties, such as the area, it is not clear how the intimate
relationship between geometrical and mechanical properties manifests.

Here we show that the connection between energy and shape of a particle
induces a novel shear weakening crossover and unusual memory properties in
various models that differ in the shape dependence of the energy of the
particles. These models exhibit distinct rheological responses at small
strains. However, all of them weaken at large strains and their shear stress
becomes negligible. During this weakening crossover, particles acquire model-
independent anisotropic shapes, orient analogously in space, and develop
short-range tetratic order [7, 8, 9, 10]. In the strain weakened regime, all
models exhibit unconventional memory properties in that under strain reversal,
the microscopic dynamics evolves irreversibly, while the contact network
changes reversibly. These results demonstrate that the link between energy and
geometry induces, in various models, a novel rheological response and
unconventional memory properties.

Figure 1: (a) Dependence of the mechanical energy (red) and shear stress
(blue), (b) the transverse transverse (red) and longitudinal (blue)
diffusivity, and (c) the eccentricity (red) and nematic order parameter (blue)
on the shear strain, for all the models. Gray and yellow shaded regions
identify the strain weakening crossover, and the strain weakened regime.

##  II Models and Methods

We consider models suitable to describe systems of deformable particles
compressed to high density that tessellate the space in two dimensions. The
centres of the particles ${\bf r}_{i}$ are the configurational degrees of
freedom. We associate a cell to each centre through a Voronoi tessellation,
and assume the cells to correspond to the particles’ shapes. We numerically
investigate four models describing systems of deformable particles at density
$\rho=1$ at which the average particle area is $1$. The energy depends on how
area $a_{i}$, perimeter, $p_{i}$ and center ${\bf r}_{i}$ of each particle
respectively differ from the average area, a dimensionless perimeter $p_{0}$,
and the centre of mass ${\bf r}_{i}^{(c)}$,

$U=\sum_{i=1}^{N}k_{a}(a_{i}-1)^{2}+k_{p}(p_{i}-p_{0})^{2}+\epsilon({\bf
r}_{i}-{\bf r}_{i}^{(c)})^{2},$ (1)

with $k_{a}$, $k_{p}$ and $\epsilon$ parameters. We study (AP) the
area+perimeter model, corresponding to the Voronoi model for epithelial cell
tissue [11, 12, 4, 13, 14], $k_{a}=1$, $k_{p}=1$, $\epsilon=0$; (C) the
centroid or Voronoi liquid model [15, 16, 17], $k_{a}=0$, $k_{p}=0$,
$\epsilon=1$; (A) the area model, $k_{a}=1$, $k_{p}=0$, $\epsilon=0$; (CA) the
centroid+area model, $k_{a}=1$, $k_{p}=0$, $\epsilon=1$. We present results
for $p_{0}=3.4$, but have verified that results only qualitatively depend on
this choice as long as $p_{0}\leq 3.81$, where the system exhibits a solid
response at small strains [11, 4]. Neglecting the Voronoi constraint on the
shape of the particles, models AP and CA impose two constraints on each
particles. Hence, the number of constraints equals the number of degrees of
freedom, and these models are marginal. Conversely, models A and C are
underconstrained, and hence are expected to have no mechanical rigidity.

We investigate the response to shear deformations of a $N=2000$ particles
enclosed in a square box with length $L=\sqrt{N}$, using Lees Edwards periodic
boundary conditions [18]. We construct the Voronoi tessellations via the C++
Boost Voronoi library [19], that accurately handles degenerate vertices. All
simulations start from initial configurations generated via the minimization
of the energy of a set of points randomly placed in the simulation box.

We study the response to shear deformation in the limit of zero shear rate and
zero temperature by resorting to athermal quasi-static shear simulations
(AQS). These simulations involve the iteration of a two-step loop, a) increase
of the shear strain, $\gamma\to\gamma+\delta\gamma$; b) relaxation of the
energy of the system. For computational efficiency, the relaxation of the
system is obtained via an energy minimization algorithm rather than resorting
to a damped dynamics. In our study, we fix $d\gamma=10^{-3}$ and resort to the
conjugate gradient minimization algorithm as implemented in the GSL library
[20].

In AQS simulations, the dependence of the elastic energy on the shear strain
comprises a series of continuous elastic branches, within which the energy
increases, and is punctuated by sudden drops corresponding to failure events.
We identify these events via a standard thresholding procedure on the energy
changes in a strain increment. Due to the many-body nature of the interaction
potential, the stress is not directly accessible via the virial theorem [21].
Hence, here we evaluate the stress as $\sigma=\frac{dU}{d\gamma}$,
approximating the elastic energy $U$ with a polynomial within each elastic
branch, to reduce numerical noise.

##  III Mechanical response

We illustrate the strain dependence of the energy and stress of the different
models in figure 1(a). At small strains, $\gamma\leq 3$, the AP model exhibits
an elastic response followed by a plastic one, as common in molecular or
colloidal systems, while the A and the C model are clearly unjammed. The CA
model appears to have a small but finite stress - in agreement with the
expectation that, due to constraint counting [12], this model does not have a
true jamming transition phase.

At larger strains, the energy of all models suddenly increases. This increase
marks the beginning of a crossover regime in the rheological properties (grey
shaded area). In the AP model, during this crossover the stress drops. In the
other models, this crossover entails a shear-jamming transition [22, 23]
leading to the growth of the shear stress, followed by a stress drop. In all
cases, we assume the crossover ends when the shear stress starts fluctuating
around zero. Due to the absence of mechanical rigidity of the final states we
interpret the observed scenario as a strain weakening crossover. The strain
weakening crossover correlates with changes in the dynamical and geometric
properties of the system. Since the system is not in a steady-state, we
evaluate the dynamics through a $\gamma$-dependent diffusivity [24],

$D_{r}(\gamma)=\langle(r(\gamma+d\gamma)-r(\gamma-d\gamma))^{2}\rangle/(2d\gamma),$
(2)

where $r=x-y\gamma$ or $r=y$ is a component of the not-affine displacement
field and $d\gamma=0.2$ is the scale of observation. Figure 1(b) shows that
the diffusivity sharply drops at the beginning of the crossover in models A,C,
which undergo a shear-jamming transition, and of the AC model, which also
rigidifies. In all models, the diffusivity reaches minute values at the end of
the crossover, $\simeq 10^{-4}$. This drop in the diffusivity is the proxy of
the emergence of an almost affine response of the system to shear
deformations: only after strain increments of order $10^{4}$ the (non-affine)
separation between close particles varies of a quantity comparable to the
typical inter-particle distance.

##  IV Evolution of geometrical properties

We investigate the evolution of geometrical properties of the system
associating to each particle an eccentricity
$e=\sqrt{1-\left(\lambda_{2}/\lambda_{1}\right)^{2}}$, where
$\lambda_{1}\geq\lambda_{2}$ are the eigenvalues of the covariance matrix of
its vertices [6]. The direction of the eigenvector associated with
$\lambda_{1}$ defines the cell’s polarity, ${\bf p}={\bf u}_{1}/{\mid{\bf
u}_{1}\mid}$. The nematic order parameter is then $S=\langle
2\cos^{2}{\phi_{i}}-1\rangle$, $\phi_{i}$ being the angle between the polarity
of cell $i$ and the nematic director. The average eccentricity of the cells
grows at small strains, but then slightly drops within the strain weakening
crossover, as apparent from figure 1(c) (red). The polarities of the cells
remain roughly uncorrelated, as the nematic order parameter (blue, right axis)
does not vary significantly, if not during the initial transient of the AP
model. Hence, in the strain weakened regime cells are more elongated than in
the initial unstrained configuration, but their polarities are not aligned
over long distances.

Figure 2: Distribution of cell edge-lengths $l$ (a,b), shape index $q_{0}$
(c,d) and orientation $\phi$ (e,f), for $\gamma=0$ and $\gamma=6$
respectively.

In the strain weakened state the system rearranges easily and hence possesses
an extensive number of soft spots [25, 3, 4, 5, 6, 26, 27], structural
rearrangements induced by a small increment $s$ of the local shear stress.
Accordingly, the distribution $P(s)$ of the local stress increments needed to
trigger a rearrangement scales as $s^{\theta}$ with $\theta>0$ in the solid
phase, and $\theta=0$ in the fluid one. It is generally arduous to infer
$P(s)$ from a snapshot of the system, and hence to correlate geometrical and
mechanical properties. However, when the energy of the particles is dominated
by their perimeter, the stress needed to trigger a T1 transition leading to
the disappearance of an edge of length $l$ is $s\propto l$. Hence,
$P(s)\propto P(l)$ [5, 28]. For the parameter values we have considered, this
occurs in the AP model [2, 3, 4, 5, 6]. Indeed, figures 2(a) and (b) reveal
that the edge-length distribution of the AP model quickly vanishes as $l\to 0$
in the solid phase at small $\gamma$, while it saturates to a finite value in
the strain weakened regime, implying $\theta=0$. For the other models, $P(l)$
is not in principle informative of the degree of mechanical stability, and we
find $P(l=0)>0$ at all $\gamma$. Surprisingly, however, we find that $P(l)$,
which is model dependent at $\gamma=0$, becomes model independent in the
strain weakened regime, as in figures 2(a,b).

Figure 3: Density correlation function $g(x,y)$, for (a) AP, (b) A, (c) C and
(d) CA models, respectively. The different columns refer to $\gamma=0,3$ and
$6$, from left to right.

Other geometrical properties also appear to become model-independent in the
strain weakened regime. For instance, figures 2(c) and (d) show that the
distribution of the shape index of the particles[29, 2], as set by their
actual shape $q_{0}=\frac{p_{i}}{\sqrt{a}_{i}}$, is model dependent in the
initial configuration and becomes model independent in the asymptotic strain
weakened regime. The average $q_{0}$ value increases as the system is
strained, consistent with the increase of the average particle eccentricity.
The distribution of the orientations $\phi$ of the polarity vectors of the
cells, ${\bf p}=(\cos\phi,\sin\phi)$, which is flat in the initial
configuration (e), is model dependent at intermediate strains (not shown) and
becomes approximately model independent in the strain weakened regimes (f). In
this regime, the distribution has a well pronounced minimum at around
$\phi\simeq\pi/2$, and for all but the AP model, two asymmetric maxima at
$\phi\simeq\pi/4$ and $3\pi/4$. The results of figure 2 thus indicate that in
the strain weakened regime particles acquire similar model-independent typical
shapes, and that they are similarly oriented in space.

The above analysis suggests that all models acquire a same typical structure
in the strain weakened regime.

Figure 4: Evolution of 2d polar correlation under strain for (a) AP, (b) A,
(c) C and (d) CA model respectively. The different columns refer to
$\gamma=0,3$ and $6$, from left to right.  Figure 5: The dashed-dot line
illustrates the evolution of the persistence (a) and of energy (b) of the AP
model, as the strain slowly increases. The other lines show how these
quantities evolve as the strain is reversed, after having reached a maximum
value $\gamma_{r}$. Panel (c) illustrates the increase of the mean square
displacement as the strain decreases. At large $\gamma$ in the strain weakened
regime, the persistence is reversible while the energy is not.

To test this hypothesis, we investigate the 2d density correlation function,
$g(x,y)$, and the 2d correlation function of the cell’s polarities,

$p(x,y)=\langle\cos 2(\phi_{i}-\phi_{j})\rangle$ (3)

where the average is over all cell pairs with a distance ($x$, $y$). In
initial configurations of all models, $g(x,y)$ is radially symmetric, with
correlations more pronounced in the AP and CA models, as in the left most
column of figure 3 and figure 4. As the strain increases, the evolution of the
structure is strongly model dependent. For instance, figure 3(middle column)
clarifies that only the AP model transiently develops crystalline order over
an appreciable length scale, the underlying crystal corresponding to an
obliquely strained triangular lattice, as apparent from the location of the
first six peaks. However, strong similarities between the correlation
functions of the different models emerge in the strain weakened regime (right
column). In particular, a four pointed star region with $g\simeq-1$, is
present in all models. The star is slightly rotated with respect to the
$x$-axis, and the two almost horizontal arms are longer than the vertical
ones. Considering that square-lattice of central-force model are generally
observed to be mechanically unstable [30, 31] at zero temperature, we
speculate that the emergence of this short-ranged tetratic order [7, 8, 9], is
a signature of the existence of an extensive number of soft spots that allow
the system to flow exerting minute mechanical stress.

Figure 6: Snapshots of a small region of configurations of the AP model
reached increasing (a) and decreasing (b) the strain. Particles are colour
coded according to the value of their shape index, $q_{0}$, and the shear is
reversed at $\gamma_{r}=20$. Panel (c) illustrates the $\gamma$ dependence of
the fluctuations of $q_{0}$.

##  V Memory in the strain weakened regime

In the strain weakened regime, the system is jammed: the energy is finite and
the shear stress fluctuates about minute values (see figure 1). In addition,
the diffusivity is strongly suppressed but not fully inhibited as it occurs,
e.g., in a strain induced diffusive-absorbing phase transition of colloidal
systems [32, 33, 34, 35]. This suggests that cell-cell contacts survive large
strains and that the contact network keeps memory of its initial state. To
check this hypothesis, we define the persistence $p_{i}(\gamma)$ of particle
$i$ as the fraction of Voronoi neighbors of the $\gamma=0$ configuration which
are also Voronoi neighbors at strain $\gamma$, and study the average
persistence as the strain first increases up to $\gamma_{r}$, and then
decreases to zero. Figure 5 (a) shows that the persistence of the AP model
behaves reversibly, at large $\gamma$ in the strain weakened regime,
regardless of the value of $\gamma_{r}$. The other models behave similarly.
Since the persistence is a measure of the evolution of the topology of the
system, we understand that as the strain is reduced the system gradually
recovers previously broken contacts. This finding indicates that in the
weakened regime the system acts as a memory device. A scalar $x$ can be
recorded by straining a weakened configuration up to $\gamma\propto x$. To
read the written quantity, one should strain the system in reverse until the
topology on the starting configuration is recovered. We expect the maximum
strain the system is able to record to be of order $D_{y}^{-1}\simeq 10^{4}$
\- a too large value to be tested numerically.

While the observed phenomenology appears analogous to the memory through path
reversal [36] which characterizes, e.g., highly viscous liquids [37], there
are interesting and crucial differences. Indeed, in viscous liquids memory
emerges as the microscopic dynamics is reversible under time reversal. This
cannot occur here as, in the strain weakened regime, stress avalanches
indicative of irreversible dynamics are clearly present. To demonstrate this
point, we investigate in figure 5(b) the evolution of the energy under strain
reversal. The energy slowly increases with $\gamma$ in the strain weakened
regime, but it does not decrease when the strain is reversed. Instead, the
energy increases up to a strain $\gamma_{t}$ along the reverse path, where the
system resets itself through a gigantic failure event to a low energy random
state that is similar to the unstrained one. Similarly, we observe in figure 5
(c) the mean square displacement
$\langle(\bf{r}(\gamma)-\bf{r}(\gamma_{r}))^{2}\rangle$ also increases as
$\gamma$ decreases. This behavior does not depend on the strain $\gamma_{r}$
reached before reversing the strain increment.

To further illustrate that cell’s properties change during the reverse shear,
we provide in figure 6 partial snapshots of the system at different $\gamma$
values reached as the strain increases (a) and decreases (b), coloring each
cell according to the values of its own shape index, $q_{0}$. Configurations
with the same $\gamma$ value, reached on increasing or increasing $\gamma$,
clearly differ. These differences are quantified in panel (c), which
illustrates that configurations reached as the strain decreases have much
larger $q_{0}$ fluctuations.

##  VI Conclusion

Our results demonstrate that, at high density, the correlation between energy
and shape of the particles induces a previously unreported strain-weakening
crossover in systems of soft deformable particles. This involves a drop in the
shear stress to minute values as the saturates and particles acquire model-
independent shapes. In the weakened regime, the contact network evolves
reversibly, while the particles’ positions do not. This crossover is enabled
by particles’ deformability, that weakens the system by allowing for the
emergence of short-ranged tetratic order.

We notice that Ref. [5] investigated the AP model under shear, focusing on the
Vertex model, and did not observe any weakening crossover. In the Vertex model
the degrees of freedom are the vertices rather than the cell centers, and
topological transitions are artificially triggered according to some rule
based on a chosen threshold for the edge length. Possibly, the presence of
this rule might critically affect the response to shear, and indeed we observe
edge-lengths of order $10^{-5}$, while the smallest edged reported in Ref. [5]
is of order $10^{-2}$. Alternatively, it is possible that the observed
phenomenology stems from the Voronoi constraint shared by all our models, i.e.
from the fact that the particles’ shape are identified with the Voronoi cells.
It is not obvious that this should be the case, as the Voronoi constraint
should stiffen the system, rather than promoting its weakening [12].

So far, our reported strain weakening crossover has not been observed in
experiments. For instance, when subjected to an imposed strain, foams do not
weaken homogeneously as we observe here, but rather develop a shear band [38,
39]. There could exist systems exhibiting a strain weakening crossover. The
response of epithelial tissue to shear strain deformations, for instance, has
not been experimentally investigated. Similarly, the response of high density
microgel suspensions in a 2d-geometry is unexplored. Hence, it would be
interesting to experimentally ascertain the rheological properties of these
systems, and of related systems of soft deformable particles at high density.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

We thank the Singapore Ministry of Education through the Academic Research
Fund 2019-T1-001-032 (Tier 1). The computational work for this article was
fully performed on resources of the National Supercomputing Centre, Singapore.

## References

  * Likos [2001] C. N. Likos, _Physics Reports_ , 2001, 348, 267–439.
  * Bi _et al._ [2015] D. Bi, J. H. Lopez, J. M. Schwarz and M. L. Manning, _Nat. Phys._ , 2015, 11, 1074–1079.
  * Kim _et al._ [2018] S. Kim, Y. Wang and S. Hilgenfeldt, _Physical Review Letters_ , 2018, 120, 248001.
  * Li and Ciamarra [2018] Y.-W. Li and M. P. Ciamarra, _Phys. Rev. Materials_ , 2018, 2, 045602\.
  * Popović _et al._ [2021] M. Popović, V. Druelle, N. A. Dye, F. Jülicher and M. Wyart, _New Journal of Physics_ , 2021, 23, 033004.
  * Li _et al._ [2021] Y.-W. Li, L. L. Y. Wei, M. Paoluzzi and M. P. Ciamarra, _Physical Review E_ , 2021, 103, 022607.
  * Donev _et al._ [2006] A. Donev, J. Burton, F. H. Stillinger and S. Torquato, _Physical Review B_ , 2006, 73, 054109.
  * Hou _et al._ [2020] Z. Hou, Y. Zong, Z. Sun, F. Ye, T. G. Mason and K. Zhao, _Nature communications_ , 2020, 11, 1–11.
  * Geng _et al._ [2009] J. Geng, J. V. Selinger _et al._ , _Physical Review E_ , 2009, 80, 011707.
  * Li _et al._ [2019] Y.-W. Li, Z.-Q. Li, Z.-L. Hou, T. G. Mason, K. Zhao, Z.-Y. Sun and M. P. Ciamarra, _Physical Review Materials_ , 2019, 3, 125603.
  * Bi _et al._ [2016] D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, _Phys. Rev. X_ , 2016, 6, 021011.
  * Sussman and Merkel [2018] D. M. Sussman and M. Merkel, _Soft matter_ , 2018, 14, 3397–3403.
  * Pasupalak _et al._ [2020] A. Pasupalak, L. Yan-Wei, R. Ni and M. P. Ciamarra, _Soft Matter_ , 2020, 16, 3914–3920.
  * Zheng _et al._ [2020] Y. Zheng, Y. W. Li and M. P. Ciamarra, _Soft Matter_ , 2020, 16, 5942–5950.
  * Hain _et al._ [2020] T. M. Hain, M. A. Klatt and G. E. Schröder-Turk, _The Journal of chemical physics_ , 2020, 153, 234505.
  * Ruscher _et al._ [2017] C. Ruscher, J. Baschnagel and J. Farago, _Phys. Rev. E_ , 2017, 97, 032132.
  * Klatt _et al._ [2019] M. A. Klatt, J. Lovrić, D. Chen, S. C. Kapfer, F. M. Schaller, P. W. Schönhöfer, B. S. Gardiner, A. S. Smith, G. E. Schröder-Turk and S. Torquato, _Nature Communications_ , 2019, 10, 1–9.
  * Lees and Edwards [1972] A. Lees and S. Edwards, _Journal of Physics C: Solid State Physics_ , 1972, 5, 1921.
  * Schäling [2011] B. Schäling, _The boost C++ libraries_ , Boris Schäling, 2011.
  * Galassi [2018] M. e. a. Galassi, _GNU Scientific Library Reference Manual_ , 2018.
  * Zimmerman _et al._ [2004] J. A. Zimmerman, E. B. WebbIII, J. Hoyt, R. E. Jones, P. Klein and D. J. Bammann, _Modelling and simulation in materials science and engineering_ , 2004, 12, S319.
  * Bi _et al._ [2011] D. Bi, J. Zhang, B. Chakraborty and R. P. Behringer, _Nature_ , 2011, 480, 355–358.
  * Pica Ciamarra _et al._ [2011] M. Pica Ciamarra, R. Pastore, M. Nicodemi and A. Coniglio, _Physical Review E_ , 2011, 84, 041308.
  * Chattoraj _et al._ [2011] J. Chattoraj, C. Caroli and A. Lemaître, _Phys. Rev. E_ , 2011, 84, 011501.
  * Lemaître and Caroli [2007] A. Lemaître and C. Caroli, _arXiv preprint arXiv:0705.3122_ , 2007\.
  * Karmakar _et al._ [2010] S. Karmakar, E. Lerner and I. Procaccia, _Physical Review E_ , 2010, 82, 055103.
  * Lin _et al._ [2014] J. Lin, A. Saade, E. Lerner, A. Rosso and M. Wyart, _EPL (Europhysics Letters)_ , 2014, 105, 26003.
  * Krajnc _et al._ [2018] M. Krajnc, S. Dasgupta, P. Ziherl and J. Prost, _Physical Review E_ , 2018, 98, 022409.
  * Park _et al._ [2015] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim, B. Gweon _et al._ , _Nature materials_ , 2015, 14, 1040–1048.
  * Żabińska [1991] K. Żabińska, _Phys. Rev. B_ , 1991, 43, 3450–3459.
  * Mao _et al._ [2015] X. Mao, A. Souslov, C. I. Mendoza and T. C. Lubensky, _Nature Communications 2015 6:1_ , 2015, 6, 1–8.
  * Gadala‐Maria and Acrivos [1980] F. Gadala‐Maria and A. Acrivos, _Journal of Rheology_ , 1980, 24, 799–814.
  * Tsai _et al._ [2003] J.-C. Tsai, G. A. Voth and J. P. Gollub, _Phys. Rev. Lett._ , 2003, 91, 064301.
  * Grebenkov _et al._ [2008] D. S. Grebenkov, M. Pica Ciamarra, M. Nicodemi and A. Coniglio, _Physical Review Letters_ , 2008, 100, 078001.
  * Zheng _et al._ [2021] Y. Zheng, A. D. S. Parmar and M. P. Ciamarra, _Physical Review Letters_ , 2021, 126, 118003.
  * Keim _et al._ [2019] N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry and S. R. Nagel, _Reviews of Modern Physics_ , 2019, 91, 035002.
  * Taylor [1986] G. I. Taylor, _Low Reynolds Number Flows_ , 1986.
  * Debrégeas _et al._ [2001] G. Debrégeas, H. Tabuteau and J.-M. D. Meglio, _Phys. Rev. Lett._ , 2001, 87, 178305.
  * Kabla and Debrégeas [2003] A. Kabla and G. Debrégeas, _Phys. Rev. Lett._ , 2003, 90, 258303\.