# Data-driven discovery of active nematic hydrodynamics

Chaitanya Joshi<EMAIL_ADDRESS>Department of Physics, Brandeis
University, Waltham, Massachusetts 02453, USA  Department of Physics and
Astronomy, Tufts University, Medford, Massachusetts 02155, USA  Sattvic Ray
Department of Physics, University of California at Santa Barbara, Santa
Barbara, California 93106, USA  Linnea M. Lemma  Department of Physics,
Brandeis University, Waltham, Massachusetts 02453, USA  Department of Physics,
University of California at Santa Barbara, Santa Barbara, California 93106,
USA  Minu Varghese  Department of Physics, Brandeis University, Waltham,
Massachusetts 02453, USA  Department of Physics, University of Michigan, Ann
Arbor, Michigan 48109 USA  Graham Sharp  Department of Physics, University of
California at Santa Barbara, Santa Barbara, California 93106, USA  Zvonimir
Dogic<EMAIL_ADDRESS>Department of Physics, University of California
at Santa Barbara, Santa Barbara, California 93106, USA  Aparna Baskaran
<EMAIL_ADDRESS>Department of Physics, Brandeis University, Waltham,
Massachusetts 02453, USA  Michael F. Hagan<EMAIL_ADDRESS>Department of
Physics, Brandeis University, Waltham, Massachusetts 02453, USA

###### Abstract

Two-dimensional active nematics are often modeled using phenomenological
continuum theories that describe the dynamics of the nematic director and
fluid velocity through partial differential equations (PDEs). While these
models provide a statistically accurate description of the experiments, the
identification of the relevant terms in the PDEs and their parameters is
usually indirect. Here, we adapt a recently developed method to automatically
identify optimal continuum models for active nematics directly from the
spatio-temporal director and velocity data, via sparse fitting of the coarse-
grained fields onto generic low order PDEs. We test the method extensively on
computational models, and then apply it to data from experiments on
microtubule-based active nematics. Thereby, we identify the optimal models for
microtubule-based active nematics, along with the relevant phenomenological
parameters. We find that the dynamics of the orientation field are largely
governed by its coupling to the underlying flow, with free-energy gradients
playing a negligible role. Furthermore, by fitting the flow equation to
experimental data, we estimate a key parameter quantifying the ‘activity’ of
the nematic.

###### pacs:

Valid PACS appear here

Active nematics demonstrate how energy-consuming motile constituents can self-
organize into diverse non-equilibrium dynamical states Marchetti  _et al._
(2013); Ramaswamy (2010); Toner  _et al._ (2005). They offer a versatile
platform to both advance our fundamental understanding of non-equilibrium
physics and develop materials with properties that are thermodynamically
forbidden in equilibrium. These twin goals require theoretical models that can
reveal the mechanism underlying the emergent dynamics, and guide rational
design to elicit desired spatio-temporal dynamics. In this paper, we combine
data-driven model discovery with experiments and computational modeling to
identify the most parsimonious model for a specific experimental realization
of active nematics. Using the discovered model, we identify the relationship
between key theoretical parameters, such as the magnitude of activity, and
experimental control variables. The described methods can be applied to
diverse realizations of active nematics ranging from shaken rods to motile
cells Narayan  _et al._ (2007); Wensink  _et al._ (2012); Zhou  _et al._
(2014); Duclos  _et al._ (2017); Kawaguchi  _et al._ (2017); Kumar  _et al._
(2018), as well as other realizations of active matter.

Our target is to describe microtubule-based active nematics. Being
reconstituted from tunable and well-characterized components, they afford a
unique opportunity to develop continuum theory models and connect these to
microscopic dynamics Sanchez  _et al._ (2012); DeCamp  _et al._ (2015);
Doostmohammadi  _et al._ (2018). Hydrodynamic theories based on the expansion
of symmetry-allowed terms have provided insight into dynamics of active
nematics in general, and the microtubule-based system specifically. For
example, such models have been used to describe defect dynamics Giomi  _et
al._ (2013, 2014); Doostmohammadi  _et al._ (2017); Oza and Dunkel (2016);
Cortese  _et al._ (2018); Shendruk  _et al._ (2018), induced flows in the
suspending fluid Thampi  _et al._ (2014); Giomi (2015); Lemma  _et al._
(2019), and how confinement in planar Shendruk  _et al._ (2017); Norton  _et
al._ (2018); Gao  _et al._ (2017) and curved geometries Zhang  _et al._
(2016); Ellis  _et al._ (2017); Alaimo  _et al._ (2017) controls defect
proliferation and dynamics. These efforts employed a range of hydrodynamic
models, but each assumed different symmetry-allowed terms. Thus, the field
lacks a standard model and understanding of magnitudes and sources of error in
existing models.

Improving hydrodynamic theories is a formidable task. Accounting for higher-
order effects involves an intractable number of additional symmetry-allowed
terms, and the phenomenological coefficient associated with each term must be
estimated from experimental data. To overcome these limitations, we identify
the most parsimonious model that captures the experimentally observed
dynamics. We build on a previous framework Brunton  _et al._ (2016); Rudy  _et
al._ (2017) that has recently been applied to particle-based simulations of
active matter Maddu  _et al._ (2022) and computational and experimental data
of overdamped polar particles Supekar  _et al._ (2021). We employ extensive
PolScope and fluorescence measurements of microtubule alignment and PIV
measurements of velocities, to identify equations governing both the
orientational dynamics and the activity-driven flows coupled to long-range
hydrodynamic interactions in the fluid momentum equation.

Figure 1: Benchmarking model discovery framework on continuum simulations
data. _Optimality curves_ : $1-R^{2}$ as a function of number of non-zero
terms in the model for (a) the $Q_{xx}$ component of Eq. (1), (b) the Stokes
equation and (c) the unsteady Stokes equation. (a): Starting with the full
library, the least important term is eliminated one by one to get a hierarchy
of models. The $R^{2}$ deteriorates when terms used in the simulation begin to
drop. (b), (c): same as (a) for the flow equation, Eq. (4), excluding the
pressure.

Our framework is built on the Sparse Identification of Non-linear Dynamics
(SINDy) algorithm Brunton  _et al._ (2016) and its extension to PDEs (PDE-
FIND) Rudy  _et al._ (2017). We adapt some key improvements of this method to
suit the microtubule-based active nematics system Alves and Fiuza (2020);
Reinbold and Grigoriev (2019); Reinbold  _et al._ (2020, 2021). With the
orientation and velocity data available from the experiments, we seek models
describing the active nematic as a single 2D fluid with nematic symmetry Beris
and Edwards (1994); Doostmohammadi  _et al._ (2018). Hence, there are two
fields: the nematic tensor order parameter
$\mathbf{Q}=s[\mathbf{n}\otimes\mathbf{n}-(1/2)\mathbf{I}]$ and a flow field
$\mathbf{u}$, with $\mathbf{n}$ as the local orientation unit vector and $s$
the scalar order parameter. $\mathbf{Q}$ is symmetric and trace-less by
definition. We assume a constant density and an incompressible fluid; the
latter is validated by numerical measurements of the divergence of the
velocity field Lemma  _et al._ (2021). Our model then consists of 4
independent scalar fields: $Q_{xx}$, $Q_{xy}$, $u_{x}$, $u_{y}$, and a latent
variable $P$ (pressure).

We begin by postulating the generalized form of the underlying model. The
dynamics of the Q-tensor takes the form common throughout all continuum
theories of active nematics:

$\partial_{t}Q_{ij}=\sum_{k}a^{k}_{ij}F_{k}(\mathbf{Q},\mathbf{u},\nabla\mathbf{Q},\nabla\mathbf{u},\ldots)$
(1)

where the $F_{k}$’s are combinations (potentially non-linear) of $\mathbf{Q}$,
$\mathbf{u}$, and their spatial derivatives up to a maximum order, and
$a^{k}_{ij}$’s are the corresponding phenomenological coefficients. For
instance, in 2D, a well-known model for the Q-equation is Doostmohammadi  _et
al._ (2018):

$\partial_{t}\mathbf{Q}+\mathbf{u}\cdot\nabla\mathbf{Q}-\mathbf{S}=\Gamma_{r}\mathbf{H}$
(2)

where
$\mathbf{S}=-(\mathbf{\Omega}\cdot\mathbf{Q}-\mathbf{Q}\cdot\mathbf{\Omega})+\lambda\mathbf{E}-2\lambda\mathbf{Q}(\mathbf{Q}\colon\nabla\mathbf{u})$
is the co-rotation term and
$\mathbf{H}=a_{2}\mathbf{Q}+a_{4}\Tr(\mathbf{Q}^{2})\mathbf{Q}+K\nabla^{2}\mathbf{Q}$
is the negative gradient of the liquid crystal free energy. Here,
$E_{ij}=(\partial_{i}u_{j}+\partial_{j}u_{i})/2$ and
$\Omega_{ij}=(\partial_{i}u_{j}-\partial_{j}u_{i})/2$ are the strain rate and
vorticity tensors respectively, $\lambda$ is the flow alignment parameter,
$\Gamma_{r}$ is the rotational diffusion coefficient, $K$ is the bending
modulus, and $a_{2}>0$, $a_{4}<0$ are phenomenological coefficients
corresponding to the isotropic-nematic transition (see Appendix Section A for
an extended discussion of the theory). We build a vast library of the terms
$F_{k}$ that can capture models well beyond Eq. (2). Further, we make no
physics-based simplifying assumptions, e.g. translational, rotational, and
Galilean invariance, for the alignment equation (Eq. 1). Hence, the discovery
of a model which satisfies these conditions is a test of the algorithm. This
results in a large number $n=246$ of candidate terms (Appendix Section B).

For the flow equation, the usual form assumed for model-discovery is Navier-
Stokes-like, with the time-derivative on the left side and rest of the terms
on the right side Rudy  _et al._ (2017); Reinbold and Grigoriev (2019);
Reinbold  _et al._ (2020, 2021). However, this approach is not appropriate for
our system. Since it is in the low Reynolds number regime Giomi (2015);
Doostmohammadi  _et al._ (2017), the significance of the time-derivative term
itself needs to be investigated. Indeed, the active nematic flow has been
modeled using pure Stokes Varghese  _et al._ (2020); Zhou  _et al._ (2021);
Duclos  _et al._ (2020), unsteady Stokes Giomi  _et al._ (2012); Giomi (2015)
as well as full Navier Stokes Giomi  _et al._ (2013); Giomi and Desimone
(2014); Doostmohammadi  _et al._ (2018); Shendruk  _et al._ (2017);
Chandragiri  _et al._ (2019, 2020); Thampi  _et al._ (2013, 2014, 2015, 2016)
formulations. While these approaches have been compared numerically Koch and
Wilczek (2021), there has yet to be a definitive indication of the
contributions of the inertial terms for this system. Since the viscous forcing
is guaranteed to exist in this regime, we assume a form

$\nabla^{2}\mathbf{u}=c_{0}\partial_{t}\mathbf{u}+\sum_{i}c_{i}\mathbf{H}_{i}(\mathbf{Q},\mathbf{u},\nabla\mathbf{Q},\nabla\mathbf{u},\ldots)$
(3)

with $\nabla\cdot\mathbf{u}=0$, and the time-derivative on the _right hand
side_ so that its contribution can be evaluated. For instance, the lowest
order symmetry-allowed ‘active stress’ in the flow equation is the well-known
$-\alpha\mathbf{Q}$, with $\alpha>0$ being the extensile ‘activity’
Doostmohammadi  _et al._ (2018); Aditi Simha and Ramaswamy (2002). In our
model form, this gives a general flow equation:

$\nabla^{2}\mathbf{u}=c_{0}\partial_{t}\mathbf{u}+c_{1}\mathbf{u}\cdot\nabla\mathbf{u}+c_{2}\nabla
P+c_{3}\nabla\cdot\mathbf{Q}+\ldots$

which can be captured by our method. In this form, the coefficient $c_{3}$
corresponds to the ratio of the activity to the viscosity, $\alpha/\eta$. In
practice, we use the stream function formulation of the flow equation to
eliminate the pressure, which is not an observable in the experiments Rudy
_et al._ (2017); Reinbold and Grigoriev (2019).

We perform model discovery from the data as follows. Setting $N_{x}$, $N_{y}$,
and $N_{t}$ as the number of measurements in the two spatial dimensions and
time respectively, we randomly select $m$ of the total $N_{x}N_{y}N_{t}$
space-time points. At each selected space-time point we evaluate a linear
system, e.g. for the $Q_{xx}$ equation, $(\partial_{t}Q_{xx})_{m\times
1}=F_{m\times n}\cdot\vec{a}_{n\times 1}$.

The derivatives must be computed numerically, which amplifies noise in the
data. To mitigate noise, we use two different approaches. In the _integral
formulation_ , for each of the $m$ selected space-time points and $n$ terms,
we compute a local average in space and time in a small window (e.g.
$\mathrm{5}\mathrm{x}\mathrm{5}\mathrm{x}\mathrm{5}$ pixels) Alves and Fiuza
(2020). This approach is effective for model discovery, but leads to
inaccurate parameter estimates for the flow equation because the noise
mitigation is not sufficient to counter noise amplification by the high-order
derivatives in that equation.

To obtain more powerful noise mitigation at the cost of additional analytical
effort, we adapt a _weak formulation_ of the PDE regression problem Reinbold
and Grigoriev (2019); Reinbold  _et al._ (2020, 2021). Briefly, we fit the
data to the weak form of Eq. (3):

$\int_{\Omega_{k}}\mathbf{w}\cdot\left[\nabla^{2}\mathbf{u}=c_{0}\partial_{t}\mathbf{u}+c_{1}\mathbf{u}\cdot\nabla\mathbf{u}+\ldots\right]$
(4)

By choosing an appropriate test function $\mathbf{w}$, we can move the
derivatives from the noisy experimental data to the exact test functions, and
also integrate out latent variables using integration by parts (see Appendix
Section B for details, the terms include in the library are available in Table
1).

Next, we seek optimal fits to these equations with the minimum number of non-
zero terms, thus yielding an interpretable model that accurately describes the
data but avoids overfitting. To this end, we perform optimization using Ridge
regression (least-squares gives similar results), and then eliminate the least
important terms one by one to obtain a hierarchy of models Alves and Fiuza
(2020). The optimal model is then determined by comparing the trade-off
between the accuracy ($R^{2}$ of the fit) and the complexity (number of non-
zero terms in the model) Alves and Fiuza (2020).

We first test our framework by recovering governing equations from data
generated by simulating active nematic hydrodynamics. To generate the data, we
consider two qualitatively different models for flow: one is purely Stokesian
with substrate friction, and the other is unsteady Stokes flow, which omits
the convective term but keeps $\partial_{t}\mathbf{u}$. Although both models
capture the primary physics of active nematics, they test the algorithm with
different types and numbers of terms Doostmohammadi  _et al._ (2018); Duclos
_et al._ (2020); Zhou  _et al._ (2021); Varghese  _et al._ (2020); Giomi  _et
al._ (2012); Giomi (2015). We solve these equations numerically and extract
the steady state Q-tensor and velocity values (Appendix Section C and Table
2). To test noise mitigation, we add 5% synthetic noise to the simulation
data, and apply the integral formulation to the alignment equation and the
weak formulation to the flow equation.

While eliminating terms one by one, we obtain the $R^{2}$ value at each model
hierarchy. We plot the _optimality curve_ as the logarithm of the Fraction of
Variance Unexplained (FVU), or ($1-R^{2}$) as a function of the number of non-
zero terms left in the model. We define the optimal number of terms $n^{*}$ as
the $n$-value at which the second derivative of the curve is highest,
indicating the largest drop in $\log(\mathrm{FVU})$ (Fig. 1). The equation
used for the alignment tensor Eq. (17), in particular for $Q_{xx}$, has 10
terms (Appendix Section C), while the Stokes and Unsteady Stokes equations
both have 2 terms. In all these cases, the framework returns the correct
equation with very small errors in the identified coefficients (Fig. 1 and
Appendix Fig. 4). Thus, we estimate important phenomenological parameters
directly from the data, including the the activity level $\alpha$, bending
modulus $K$, flow alignment coupling $\lambda$, and bulk free energy
coefficients $a_{2}$ and $a_{4}$. We benchmark the noise mitigation strategies
by performing the same protocol with increasing window sizes (Appendix Fig. 5)
and varying levels of synthetic noise (Appendix Fig. 5). The results show that
the integral formulation with a small window size of a few pixels is suitable
for the alignment equation, while the weak formulation with a large window
size, comparable to the system size, is better for the flow equation.

Figure 2: Discovering active nematic hydrodynamics from experimental data. (a)
Optimality curves (($1-R^{2}$) vs number of non-zero terms) for the $Q_{xx}$
equation from experiments with various ATP concentrations. The shoulder at
$n\sim 6$ contains flow coupling terms — convection, rotation and flow
alignment. Higher order flow alignment terms show up at higher $n$. The blue,
orange, green and red colors correspond to
$25\text{\,}\mathrm{\SIUnitSymbolMicro M}$,
$50\text{\,}\mathrm{\SIUnitSymbolMicro M}$,
$100\text{\,}\mathrm{\SIUnitSymbolMicro M}$ and
$500\text{\,}\mathrm{\SIUnitSymbolMicro M}$ ATP respectively, in (a-c). (b)
Values of the coefficients of key flow-coupling terms appearing in the optimal
models for various ATP concentrations. (c) Optimality curves for the weak form
flow equation. The activity term $\nabla\cdot\mathbf{Q}$ dominates, suggesting
Stokesian dynamics. (d) The fitted coefficient of the activity term,
$\alpha/\eta$, as a function of the ATP concentration (green circles). This
quantity compares well with the inverse of the velocity correlation time (blue
triangles), suggesting that $\alpha/\eta$ corresponds to a relevant timescale
in the system.

Next, we perform model discovery on orientation and velocity data extracted
from fluorescence images of microtubule-based active nematic (Appendix Section
D). We varied the ATP concentration, which determines the motor stepping speed
and thus determines the structure and dynamics of active nematics. Figs. 2a,c
show the optimality curves for the alignment and flow equations respectively
at various ATP concentrations. In most cases, we obtain the optimal model:

$\displaystyle\partial_{t}\mathbf{Q}\sim\ $
$\displaystyle-\mathbf{u}\cdot\nabla\mathbf{Q}-(\mathbf{\Omega}\cdot\mathbf{Q}-\mathbf{Q}\cdot\mathbf{\Omega})$
$\displaystyle+\mathbf{E}-2(\mathbf{Q}\colon\nabla\mathbf{u})\mathbf{Q}$
$\displaystyle\eta\nabla^{2}\mathbf{u}\sim\ $
$\displaystyle+\alpha\nabla\cdot\mathbf{Q}+\nabla P$ (5)

As tests of consistency, the alignment equation recovers Galilean invariance
from the data: the convective and co-rotational derivatives have coefficients
of $\sim 1$ (Fig. 2b). Furthermore, the flow alignment parameter, $\lambda\sim
1$ (Fig. 2b), is consistent with the theoretical result for the high aspect
ratio $a\gg 1$ of the microtubules, $\lambda=(a^{2}-1)/(a^{2}+1)\to 1$ Maitra
_et al._ (2018) 111A higher accuracy result for the Q-tensor equation was
obtained using a dataset with a small field of view with a few defects per
frame (see Appendix Section E). That dataset however doesn’t give a good
$R^{2}$ for the flow due to the limitation on the window size.. Interestingly,
the terms arising from the bulk liquid crystal free energy (see Eq. (2)) are
absent, a finding that supports a previous model Thampi  _et al._ (2015).
Elastic free energy terms such as $\nabla^{2}\mathbf{Q}$ are also absent.
These results indicate that the alignment dynamics of microtubule-based active
nematics are dominated by flow coupling. In comparison, contributions from the
free energy dissipation to the dynamics are negligible. The optimality curves
for the flow equation are almost flat (Fig. 2c), showing that the activity
term $\alpha/\eta\nabla\cdot\mathbf{Q}$ alone dominates the dynamics. The
inertial terms are absent (not appearing until $n^{*}\sim 5$), indicating that
the Stokes flow approximation accurately describes the experimental active
nematic. Finally, the absence of the substrate friction term indicates that
the flows are largely unscreened.

Figure 3: Comparison of the simulated optimal model with the experimental
data. (a) Inverse of the lifetime of $+1/2$ defects plotted as a function of
the value of $\alpha/\eta$ obtained from the optimal model (green circles) and
used in the simulations (blue squares). The height of the errorbars is twice
the standard error of mean. (b) Defect density plotted against the value of
$\alpha/\eta$ obtained from the optimal model (green circles) and used in the
simulations (blue squares). The density in the simulations is scaled by a
constant. The height of the errorbars is twice the standard deviation.

The terms in the optimal model for $\mathbf{Q}$ have dimensionless
coefficients. Thus, we cannot extract length scales or time scales from the
Q-tensor equation, or provide quantitative estimates for dimensional
experimental parameters characterizing $\mathbf{Q}$ dynamics 222We can force
the framework to estimate dimensional parameters by constraining the
regression procedure to include specific terms, while performing sparse
regression on the remaining terms. For example, by forcing a term
$K\nabla^{2}Q_{xx}$ (see Eq. (2)), we obtain a value for the elastic modulus
of $K\sim$ $1\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{2}\text{\,}{\mathrm{s}}^{-1}$. However, because this term has a
negligible contribution to the dynamics of $\mathbf{Q}$, the quantitative
accuracy of this estimate may be limited., which in the context of the
discovered optimal model, is an emergent consequence of the dynamics. The fit
to the flow equation provides a direct estimate of the scaled activity
parameter $\alpha/\eta$, an intrinsic time scale Giomi  _et al._ (2011), as a
function of the ATP concentration (Fig. 2d) 333Independent measurements of the
viscosity such as reported in Guillamat  _et al._ (2016) can be then used to
estimate $\alpha$, the strength of the active force. Determining the
relationship between activity and experimental control parameters has been a
significant challenge Lemma  _et al._ (2019). To test this estimate against
another observable, we compare this time scale to the velocity autocorrelation
time $\tau$, defined as $\bar{C}_{v}(\tau)=1/e$, with the autocorrelation
function
$\bar{C}_{v}(t)=\expectationvalue{\expectationvalue{\mathbf{u}(\mathbf{r},t^{\prime}+t)\cdot\mathbf{u}(\mathbf{r},t^{\prime})}_{t^{\prime}}/\expectationvalue{\mathbf{u}(\mathbf{r},t^{\prime})\cdot\mathbf{u}(\mathbf{r},t^{\prime})}_{t^{\prime}}}_{\mathbf{r}}$.
These observables closely agree over the whole range of ATP concentrations
(Fig. 2d).

Finally, we test the optimal model by performing simulations of Eq. (5) 444For
numerical stability, we include the $K\nabla^{2}\mathbf{Q}$ term in the
$\mathbf{Q}$ equation with $K=1$ and comparing the mean defect lifetime and
defect density as a function of $\alpha/\eta$ (Fig. 3). Remarkably, the defect
lifetimes for experiments and simulations align well without any pre-factors
(Fig.3a), indicating a scale-free mapping between the time units. The defect
densities from experiment and simulation also match, within a constant scaling
factor that cannot be directly estimated since we lack a length scale (Fig.3b)
555We re-scale the defect density so that the simulation value for
$\alpha/\eta=0.3$ falls on the linear fit of the experimental values. .

In summary, we have applied a data-driven PDE discovery method to identify
equations governing both the orientational dynamics and the activity-driven
flows of microtubule-based active liquid crystals. The optimal model is
surprisingly minimal. It demonstrates that: (1) flow coupling dominates the
orientational dynamics, and (2) the lowest-order active stress, proportional
to the local orientational order, together with the vanishing Reynolds number
limit describe the flow. The identified equations enable mapping between key
model parameters and experimental control variables, including the elusive
relationship between the magnitude of activity and ATP concentration. Through
comparison of several noise reduction approaches on the experimental data, as
well as extensive benchmarking against simulated data, we have identified an
approach to model discovery which is highly robust against experimental noise.
We have delineated associated requirements for experimental data acquisition;
i.e., requirements on the amount of data as well as its temporal and spatial
resolution. The approach is general, and can be applied to a wide variety of
active matter systems as well as other experimental and computational
observations of dynamics.

###### Acknowledgements.

We acknowledge support from the National Science Foundation (NSF) DMR-1855914,
OAC-2003820 and the Brandeis Center for Bioinspired Soft Materials, an NSF
MRSEC (DMR-2011486), as well as computing resources through NSF XSEDE
allocation TG-MCB090163 (Stampede and Comet) and the Brandeis HPCC which is
partially supported by the NSF through DMR-MRSEC 2011486 and OAC-1920147. We
thank Link Morgan for providing early experimental data for testing, Saaransh
Singhal for providing the simulation data for the unsteady Stokes equation,
Peter J. Foster for providing feedback on the manuscript, and Michael M.
Norton and Seth Fraden for valuable discussions.

## Appendix A: Continuum Theory of Active Nematics

The theory of active nematic suspensions builds on equilibrium nematic
hydrodynamics de and Prost (1993); Beris and Edwards (1994) and extends it to
include non-equilibrium ‘active stresses’ in the fluid Ramaswamy (2010);
Marchetti  _et al._ (2013); Giomi  _et al._ (2012, 2014); Thampi  _et al._
(2014); Gao  _et al._ (2015); Giomi (2015); Maitra  _et al._ (2018); Norton
_et al._ (2018); Hemingway  _et al._ (2016); Ngo  _et al._ (2014). The most
generic models of active nematics have considered the dynamics of not only the
orientation and the velocity, but also of the concentration Giomi  _et al._
(2012, 2013); Giomi and Desimone (2014) and/or density Thampi  _et al._
(2015). In this work, we assume uniform density and concentration everywhere
(a reasonable assumption for dense 2D bulk systems, while we do believe that
allowing the density/concentration to vary will only improve our results). The
orientational order is often defined using a tensor order parameter
$\mathbf{Q}=S[\mathbf{n}\otimes\mathbf{n}-(1/2)\mathbf{I}]$. This definition
has the $\mathbf{n}\to-\mathbf{n}$ nematic symmetry built-in, and also
provides a scalar order parameter $S(\mathbf{r},t)$ that represents the
magnitude of the orientational order at the given location. For equilibrium
systems, the free energy near the isotropic-nematic transition takes the form

$\mathcal{F}=\frac{a_{2}}{2}Q^{2}+\frac{a_{4}}{4}Q^{4}+\frac{K}{2}(\nabla
Q)^{2}$ (6)

Here, $a_{2}=(1-\rho)$ and $a_{4}=(\rho+1)/\rho^{2}$. The density $\rho$
controls the transition from the isotropic ($\rho<1$) to the nematic
($\rho>1$) state. If $\rho=\rho_{0}>1$, the minimum of the free energy is in
an ordered state, with $S_{eqm}=\sqrt{-2a_{2}/a_{4}}$. $K$ is the elastic
modulus of the nematic, assuming equal moduli for splay and bend deformations.
For the alignment tensor, the dynamical equations of motion are:

$\frac{\mathrm{D}\mathbf{Q}}{\mathrm{D}t}=\partial_{t}\mathbf{Q}+\mathbf{u}\cdot\nabla\mathbf{Q}+(\mathbf{\Omega}\cdot\mathbf{Q}-\mathbf{Q}\cdot\mathbf{\Omega})=\lambda\mathbf{E}+D_{r}\mathbf{H}$
(7)

The left hand side corresponds to the co-moving co-rotational derivative of
the Q-tensor, with
$\mathbf{\Omega}=1/2(\nabla\mathbf{u}-(\nabla\mathbf{u})^{T})$ and
$\mathbf{E}=1/2(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})$. $\lambda$ is the
flow alignment parameter, with $\lambda\mathbf{E}$ being the lowest order
contribution to the flow alignment. The next leading order term
$-2\lambda(\mathbf{Q}\colon\mathbf{E})\mathbf{Q}$ 666The second order term —
$\lambda\\{\mathbf{Q}\cdot\mathbf{E}+\mathbf{E}\cdot\mathbf{Q}-2/d\Tr(\mathbf{Q}\cdot\mathbf{E})\mathbf{I}\\}$
Doostmohammadi  _et al._ (2018); Varghese  _et al._ (2020)—vanishes
identically in 2D is also often included Doostmohammadi  _et al._ (2018).
Lastly, $\mathbf{H}=-\delta\mathcal{F}/\delta\mathbf{Q}$ and $D_{r}$ is the
coefficient of rotational diffusion. This term is the free energy
contribution, while all the others come from the coupling to the underlying
flow.

The corresponding fluid flow can be described by the incompressible Navier-
Stokes (NS) equations. These are augmented with the (passive) back-flow
$\sigma_{p}$ due to the coupling to the nematic. In addition, the lowest order
non-equilibrium nematic stress takes the form $-\alpha\mathbf{Q}$, where
$\alpha$ indicates the strength of the active forces, or “activity”. Adding
this contribution in, we get the _active_ nematic fluid equation:

$\rho(\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u})=\eta\nabla^{2}\mathbf{u}-\nabla
P+\nabla\cdot(\sigma_{p}-\alpha\mathbf{Q})$ (8)

with $\rho$ being the density and $\eta$ being the viscosity.

For thin 2D samples such as the ones considered in this paper, Eq. (8) needs
to be integrated along the height $h$ of the confinement to get a quasi-2D
equation. Integrating the viscous term gives an effective linear friction
$-\Gamma\mathbf{u}$, with $\Gamma\sim\eta/h^{2}$ Maitra  _et al._ (2018).

Since the flows have a low-Reynolds number, it is common to neglect the non-
linear convection term Giomi  _et al._ (2012); Giomi and Desimone (2014), and
often the entire inertia term Duclos  _et al._ (2020); Zhou  _et al._ (2021);
Varghese  _et al._ (2020). This is the Stokes limit. The theoretical
consequences of varying the advective inertia and substrate friction have been
the subject of numerical studies Thampi  _et al._ (2014); Koch and Wilczek
(2021).

Additionally, a higher order active stress
$\nabla(\mathbf{Q}\cdot\nabla\cdot\mathbf{Q})$ in 3D can give rise to a non-
equilibrium active force $\mathbf{Q}\cdot\nabla\cdot\mathbf{Q}$ in 2D, which
can be expected to have a similar magnitude as the primary active force
$-\alpha\nabla\cdot\mathbf{Q}$ Maitra  _et al._ (2018).

This brief sketch of active nematics models shows the wealth of information
that a data-driven method like SINDy could uncover when applied to
experimental systems. For the numerical model in this paper, we use the Stokes
limit described above. Further, we set $\sigma_{p}=0$ as it has been observed
numerically that active stresses dominate the passive stresses Thampi  _et
al._ (2016); Doostmohammadi  _et al._ (2018). We include the substrate
friction as it adds another degree of freedom that our framework has to
identify. With this, we get our dynamical equation for the fluid:

$\displaystyle\eta\nabla^{2}\mathbf{u}-\Gamma\mathbf{u}-\nabla
P-\nabla\cdot(\alpha\mathbf{Q})$ $\displaystyle=0$ (9)
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0$

The Q-tensor equation is unaltered by activity because the active stresses
come into picture through the flow coupling.

## Appendix B: Library creation and benchmarking

Each data-set is pre-processed such that the fields $Q_{xx}$, $Q_{xy}$,
$u_{x}$ and $u_{y}$ lie on the same $N_{x}\times N_{y}$ grid. With $N_{t}$
such measurements in time, we get $N_{x}\times N_{y}\times N_{t}$ values for
each of the field.

#### Integral Formulation

We begin by creating a database of terms containing the fields and their
derivatives. We compute the time and space derivatives numerically using a
central difference scheme. (Since
$\nabla\cdot\mathbf{u}=\partial_{x}u_{x}+\partial_{y}u_{y}=0$, we discard
$\partial_{y}u_{y}$ from our database.) To form the library, we then make all
multiplicative combinations of these terms with a total functional order up to
$f$ and a total gradient order up to $d$. We put further limits on the
function and gradient order of $\mathbf{u}$ appearing in the terms. Thus, we
specify two pairs, ($f_{t}$ , $d_{t}$) for the overall constraint, and
($f_{u}$, $d_{u}$) for further constraint on the terms involving velocity.
This allows us to make the computation more tractable. Motivated by
theoretical models discussed in the Appendix Section A, we use
($f_{t}=3,d_{t}=2$) and ($f_{u}=1,d_{u}=1$) for the $\mathbf{Q}$ tensor
equation (Eq. (1)) and ($f_{t}=2,d_{t}=2$) for the flow equation (Eq. (3)). In
this approach, we don’t make any simplifying assumptions about the terms : all
terms appear in an “unfolded” form, with the derivatives and inner products
expanded out.

Since our data has two spatial and one time dimensions, we have a large number
of data-points. Hence, we compute the library terms only on a sub-sample of
the data Rudy  _et al._ (2017) : we randomly select $m=5000$ points from the
$N_{x}\times N_{y}\times N_{t}$ grid and compute the local average of the
terms near the points using a small (e.g.
$\mathrm{5}\mathrm{x}\mathrm{5}\mathrm{x}\mathrm{5}$ pixels) averaging window
Alves and Fiuza (2020).

#### Weak Formulation

For the fluid flow, we assume a generalized flow equation

$\int_{\Omega_{k}}\mathbf{w}\cdot\left[\nabla^{2}\mathbf{u}=c_{0}\partial_{t}\mathbf{u}+c_{1}\mathbf{u}\cdot\nabla\mathbf{u}+\ldots\right]$

Number | Term
---|---
1 | $\partial_{t}\mathbf{u}$
2 | $\mathbf{u}\cdot\nabla\mathbf{u}$
3 | $\nabla\cdot\mathbf{Q}$
4 | $\mathbf{Q}\cdot\nabla\cdot\mathbf{Q}$
5 | $\mathbf{u}$
6 | $\mathbf{Q}\cdot\mathbf{u}$
7 | $\Tr(\mathbf{Q}^{2})\mathbf{u}$
Table 1: Terms appearing in the library for the flow equation.

Using a similar “unfolded” form for the library would generate a lot of terms,
and it is not feasible to perform the integration by parts for all of them.
Hence, we create the library by hand with a judicious choice of terms, listed
in Table 1. Here, $\mathbf{w}$ is a vector test function and the integration
domain $\Omega_{k}$ is a rectangular box of size $2H_{x}\times 2H_{y}\times
2H_{t}$, centered at $(x_{k},y_{k},t_{k})$ Reinbold  _et al._ (2020).
Following Reinbold  _et al._ (2020), we use the test function

$\mathbf{w}=\nabla\times(\psi\hat{z})$ (10)

with

$\psi=\sin(\pi\text{\@text@baccent{$t$}})(\text{\@text@baccent{$x$}}^{2}-1)^{p}(\text{\@text@baccent{$y$}}^{2}-1)^{p}$
(11)

where the underbar represents the rescaled variables
$\text{\@text@baccent{$x$}}=(x-x_{k})/H_{x}$,
$\text{\@text@baccent{$x$}}=(y-y_{k})/H_{y}$ and
$\text{\@text@baccent{$x$}}=(t-t_{k})/H_{t}$. Eq. (10) implies that
$\nabla\cdot\mathbf{w}=0$, which facilitates the elimination of pressure (see
below). The functional form of Eq. (11) guarantees that $\mathbf{w}$ and its
derivatives vanish at the domain boundaries given a sufficiently large value
of $p$. This is useful for the integration by parts, which goes as follows:

$\displaystyle u^{k}_{0}$
$\displaystyle=\int_{\Omega_{k}}\\!\mathbf{w}\cdot\nabla^{2}\mathbf{u}\
\differential\Omega$
$\displaystyle=\int_{\Omega_{k}}\\!\mathbf{u}\cdot\nabla^{2}\mathbf{w}\
\differential\Omega$ (12) $\displaystyle u^{k}_{1}$
$\displaystyle=\int_{\Omega_{k}}\\!\mathbf{w}\cdot\partial_{t}\mathbf{u}\
\differential\Omega$
$\displaystyle=-\int_{\Omega_{k}}\\!\mathbf{u}\cdot\partial_{t}\mathbf{w}\
\differential\Omega$ (13) $\displaystyle u^{k}_{2}$
$\displaystyle=\int_{\Omega_{k}}\\!\mathbf{w}\cdot(\mathbf{u}\cdot\nabla\mathbf{u})\
\differential\Omega$
$\displaystyle=-\int_{\Omega_{k}}\\!\mathbf{u}\cdot(\mathbf{u}\cdot\nabla)\mathbf{w}\
\differential\Omega$ (14)

The activity term integrates as follows,

$u^{k}_{3}=\int_{\Omega_{k}}\\!\mathbf{w}\cdot(\nabla\cdot\mathbf{Q})\
\differential\Omega=-\int_{\Omega_{k}}\\!\mathbf{Q}\colon(\nabla\mathbf{w})\
\differential\Omega$ (15)

while the pressure term integrates out to zero as follows:

$\int_{\Omega_{k}}\\!\mathbf{w}\cdot\nabla P\
\differential\Omega=-\int_{\Omega_{k}}\\!P\nabla\cdot\mathbf{w}\
\differential\Omega=0$ (16)

The derivative on $\mathbf{Q}$ in term 4 in Table 1 cannot be fully
transferred to $\mathbf{w}$ because of the non-linearity, so we have to
integrate that term directly. Similarly, terms $57$ are integrated directly as
they do not contain any derivatives.

We benchmark these methods against varying noise levels as well as varying
window sizes. To find the appropriate window size, we benchmark with a
numerical data of size $512\times 512\times 500$ with $5\%$ added noise. We
measure the $R^{2}$ value of the fit as well as the average % error in the
coefficients in the optimal model (if found correctly) as a function of the
integration window size (Fig. 4). For the integral formulation (Fig. 4 left),
we find that a small window size is sufficient to mitigate the noise. For the
experimental dataset results in Fig. 2, we use a window size of $25\times
25\times 25$ pixels, or $73\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$\times$$73\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$\times$$12\text{\,}\mathrm{s}$. In the weak formulation, larger window
sizes are needed to better sample the test function Reinbold  _et al._ (2020,
2021). We find that a window that is almost as large as the field of view in
space, and $\sim 5$ times the velocity correlation time in the time dimension
Reinbold  _et al._ (2021) works well. We perform a similar analysis with the
experimental data-sets to choose the window size. To avoid self-selection, we
use the largest of the correlation times from the data-sets to set the window
size. In Fig 2, we use a window size of $205\times 205\times 605$ pixels, or
$\sim$ $600\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$\times$$600\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$\times$$300\text{\,}\mathrm{s}$. As this window size is large, we take
$m=50$ measurements for the weak form.

We now investigate the same two parameters, namely the $R^{2}$ and the error
in coefficients, but with varying noise levels (Fig. 5). Our analyses, and
indeed our results, indicate that the integral formulation is sufficient for
the Q-tensor equation, while the weak formulation works adequately for the
flow equation.

Result:
$\mathbf{A},r2s=\texttt{HRidge}(\mathbf{F},\mathbf{\dot{X}},\lambda_{2}$)

$m$ = size($\mathbf{F}$)[2] ;

// Total number of terms

$F_{n}=\texttt{norm}(\mathbf{F})$ ;

// Column-wise norm

$\mathbf{F}_{0}=\mathbf{F}/F_{n}$ ;

// Normalize the terms for comparison

$\mathbf{A}\leftarrow$ Vector(size: $m\times m$) ;

// Store optimal model at all hierarchies

$r2s\leftarrow$ Vector(size: $m$) ;

// Store model accuracies at all hierarchies

$\mathbf{a}=\operatorname*{argmin}_{\vec{a}}\left(\lVert\dot{\mathbf{X}}-\mathbf{F}_{0}\cdot\mathbf{a}\rVert_{2}+\lambda_{2}\lVert\mathbf{a}\rVert_{2}\right)$
;

// Initial guess using Ridge

$\mathbf{A}[:,m]=\mathbf{a}/F_{n}$ ;

// Store un-normalized coefficients

$r2s[m]=\texttt{rsquared}(\mathbf{F}_{0},\mathbf{\dot{X}},\mathbf{a})$

while  _$len(\mathbf{a}) >1$_ do

$j_{0}=\operatorname*{argmin}_{j}|a_{j}|$ ;

// Find the position of the smallest coefficient

$\mathbf{a}=\mathbf{a}\setminus\\{a_{j_{0}}\\}$ ;

// Remove smallest coefficient

$\mathbf{F}_{0}=\mathbf{F}_{0}\setminus\mathbf{F}_{0}[:,j_{0}]$ ;

// Remove vector corresponding to it

$\mathbf{a}=\operatorname*{argmin}_{\vec{a}}\left(\lVert\dot{\mathbf{X}}-\mathbf{F}_{0}\cdot\mathbf{a}\rVert_{2}+\lambda_{2}\lVert\mathbf{a}\rVert_{2}\right)$;

$k=len(\mathbf{a})$;

$\mathbf{A}[:,k]=\mathbf{a}/F_{n}$;

$r2s[k]=\texttt{rsquared}(\mathbf{F}_{0},\mathbf{\dot{X}},\mathbf{a})$;

Algorithm 1 Hierarchical Ridge Regression (HRidge)

## Appendix C: Continuum Simulations

For the Q-tensor, we use the simple form (in non-dimensionalized units)

$\partial_{t}\mathbf{Q}+\nabla\cdot(\mathbf{u}\mathbf{Q})+(\mathbf{\Omega}\cdot\mathbf{Q}-\mathbf{Q}\cdot\mathbf{\Omega})=\lambda\mathbf{E}+\mathbf{H}$
(17)

with

$\displaystyle\Omega_{ij}$
$\displaystyle=\frac{1}{2}(\partial_{i}u_{j}-\partial_{j}u_{i})$
$\displaystyle E_{ij}$
$\displaystyle=\frac{1}{2}(\partial_{i}u_{j}+\partial_{j}u_{i})$
$\displaystyle H_{ij}$
$\displaystyle=(-a_{2}-a_{4}Q_{kl}Q_{lk})Q_{ij}+K\partial_{k}\partial_{k}Q_{ij}$

where $\lambda$ is the flow alignment parameter, $K$ is the bending modulus
and $a_{2}<0,a_{4}>0$ drive the system to the nematic phase. The two different
flow equations used are

$\displaystyle\eta\nabla^{2}\mathbf{u}$ $\displaystyle=\nabla
P+\Gamma\mathbf{u}+\alpha\nabla\cdot\mathbf{Q}\quad\text{(Stokes)}$ (18)
$\displaystyle\rho\partial_{t}\mathbf{u}$
$\displaystyle=\eta\nabla^{2}\mathbf{u}-\nabla
P-\alpha\nabla\cdot\mathbf{Q}\quad\text{(Unsteady Stokes)}$ (19)
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0$ (20)

Here, $\eta$ is the viscosity, $\alpha$ is the activity and $\Gamma$ is the
substrate friction. We set $\rho=1$ for the unsteady Stokes equation.

For the Stokes equation Eq. (18), we use a semi-implicit finite difference
time stepping scheme based on a convex splitting of the nematic free energy
Zhao and Wang (2016); Varghese  _et al._ (2020). To solve the Stokes equation
with incompressibility, we implement a Vanka type box smoothing algorithm on a
staggered grid Vanka (1986); Varghese  _et al._ (2020). The solution at each
time-step is found using Gauss-Seidel relaxation iterations, and the rate of
convergence to the solution is accelerated by using a multigrid method
Varghese  _et al._ (2020). The simulation codes are all in-house and are
written in C. We solve the equations in a square domain of size $200\times
200$ (in simulation units) with periodic boundary conditions. We sample the
evolution at a time step of $1$ (in simulation units) on a rectangular grid
with $\differential x\sim 0.4$. For the unsteady Stokes, the coupled equations
are numerically solved in a forward time centred space (FTCS) scheme with an
explicit time-stepped predictor-corrector method. The incompressible unsteady
Stokes equations use a modified Chorin’s projection method Hermle (2017) to
compute the velocities. The data was sampled at a time-step of $0.02$ on a
rectangular grid with $\differential x\sim 0.2$. The parameters used in both
the cases are documented in Table 2.

The equation for $Q_{xx}$ has 10 terms:

$\displaystyle\partial_{t}Q_{xx}=$ $\displaystyle-
u_{x}\partial_{x}Q_{xx}-u_{y}\partial_{y}Q_{xx}-Q_{xy}\partial_{x}u_{y}+Q_{xy}\partial_{y}u_{x}$
$\displaystyle+\lambda\partial_{x}u_{x}-a_{2}Q_{xx}-2a_{4}Q_{xx}^{3}-2a_{4}Q_{xy}^{2}Q_{xx}$
$\displaystyle+K\partial_{x}^{2}Q_{xx}+K\partial_{y}^{2}Q_{xx}$

The flow equations in the integral formulation are not used in the main text,
but they would be in the form of vorticity equations (taking the curl of Eq.
(18) and Eq. (19)):

$\displaystyle\eta\nabla^{2}\omega$
$\displaystyle=\Gamma\omega+\alpha\partial_{x}^{2}Q_{xy}-2\alpha\partial_{x}\partial_{y}Q_{xx}-\alpha\partial_{y}^{2}Q_{xy}$
$\displaystyle\eta\nabla^{2}\omega$
$\displaystyle=\partial_{t}\omega+\alpha\partial_{x}^{2}Q_{xy}-2\alpha\partial_{x}\partial_{y}Q_{xx}-\alpha\partial_{y}^{2}Q_{xy}$

However, for the weak formulation, we get two terms each, after the pressure
is eliminated as described in Appendix Section B:

$\displaystyle\eta\int_{\Omega_{k}}\mathbf{w}\cdot\nabla^{2}\mathbf{u}$
$\displaystyle=\alpha\int_{\Omega_{k}}\mathbf{w}\cdot\nabla\cdot\mathbf{Q}+\Gamma\int_{\Omega_{k}}\mathbf{w}\cdot\mathbf{u}$
$\displaystyle\eta\int_{\Omega_{k}}\mathbf{w}\cdot\nabla^{2}\mathbf{u}$
$\displaystyle=\alpha\int_{\Omega_{k}}\mathbf{w}\cdot\nabla\cdot\mathbf{Q}+\int_{\Omega_{k}}\mathbf{w}\cdot\partial_{t}\mathbf{u}$

Parameter | Stokes | Unsteady Stokes
---|---|---
$\eta$ | 1 | 1
$K$ | 1 | 1
$\alpha$ | 0.3 | 4.0
$a_{2}$ | -0.3 | -16
$a_{4}$ | 1.36 | 32
Table 2: Parameters used for the Stokes and Unsteady Stokes simulations.

## Appendix D: Experimental Methods

The active nematic samples were assembled following previously established
methods Sanchez  _et al._ (2012); DeCamp  _et al._ (2015). The active mix
consisted of microtubules, kinesin motor clusters, depleting agent, and an ATP
regeneration system. Tubulin was purified from bovine brain, labeled with
Alexa 647 dye, and polymerized in the presence of GMPCPP Castoldi and Popov
(2003); Hyman  _et al._ (1991). A truncated and biotinylated version of
Kinesin-1 (K401-BCCP-HIS) was expressed in E. Coli and purified using
immobilized metal affinity chromatography Subramanian and Gelles (2007). Motor
clusters were formed by incubating $11\text{\,}\mathrm{\SIUnitSymbolMicro L}$
of the biotinylated kinesin (0.7 mg/ml) with
$5\text{\,}\mathrm{\SIUnitSymbolMicro L}$ of streptavidin (0.35 mg/ml) on ice
for 30 min. Polyethylene glycol (35000 kDa, 1%) was used to induce microtubule
bundling. A biochemical regeneration system consisting of adenosine
triphosphate (ATP, 25$\mu$M-1.4 mM), phosphoenol pyruvate (PEP, 26 mM), and
pyruvate kinase/lactic dehydrogenase (PK/LDH) kept ATP concentration constant.
Lastly, an oxygen scavenging system consisting of glucose (0.67 mg/ml),
glucose oxidase (0.08 mg/ml), catalase (0.4 mg/ml), DTT (5.6 mM), and Trolox
(2 mM) was used to minimize sample bleaching. The components of the active mix
were combined in M2B buffer (80 mM PIPES, 1mM EGTA, 2mM MgCl2, pH 6.8).

Flow chambers were created with Parafilm sandwiched between a glass slide and
a coverslip. The glass slide was made hydrophobic with a Rain-X coating, and
the coverslip was passivated with acrylamide coating Lau  _et al._ (2009). To
assemble an active nematic, the chamber was first filled with HFE oil
containing fluoro-surfactant (0.5% w/w, RAN Biotech), followed by the active
mix. The sample was sealed with UV glue (Norland optical adhesive). The active
nematic sedimented to the oil-water interface and reached a steady state after
about an hour, and was then imaged on a spinning disk confocal microscope
using an Hamamatsu Orca-Fusion BT CMOS camera and 20$\times$ magnification.
For each sample, a sequence of 10000 images was acquired at 2 frames/sec,
except for the 25 $\mu$M ATP samples, for which 1000 frames were acquired at
0.1 frames/sec.

## Appendix E: Additional methods

### Orientation and velocity fields

The orientation and velocity fields were computed simultaneously on the
fluorescent images obtained from spinning disk confocal miscrocopy. The
orientation fields were measured using an in-house structure-tensor-based code
Rezakhaniha  _et al._ (2012); Duclos  _et al._ (2020) (written in MATLAB) on
the fluorescence images. The molecular tensor
($\mathbf{n}\mathbf{n}-1/2\mathbf{I}$) was computed from the orientation data,
and was then coarse-grained with a Gaussian smoothing filter with $\sigma=10$
pixels to obtain the Q-tensor. The velocity fields were measured using
particle-image velocimetry implemented by the MATLAB-based PIVLab software.
The fields computed by PIVLab were post-processed using a Direct Cosine
Transform - Penalized Least Squares (DCT-PLS) approach that validates the raw
data, replaces the spurious and missing vectors and does some smoothing.
Garcia (2011). The orientation fields are computed on a high resolution grid
of $1152\times 1152$ pixels, while the velocity fields are on a coarser grid
of $71\times 71$ pixels. Therefore, both fields are interpolated on an
intermediate grid of $256\times 256$ pixels.

### Defect detection and tracking

To locate the defects, we compute a map of the signed winding number
$w=1/(2\pi)\oint\nabla\theta\cdot\differential{\vec{s}}$ at every point in
space Kamien (2002); Norton  _et al._ (2018) with an integration ring of
radius of 5 pixels. The winding number is zero everywhere except at the defect
locations DeCamp  _et al._ (2015); Ellis  _et al._ (2017). To eliminate
spurious defects, we filter out regions with a non-zero winding number that
are smaller than 60 squared pixels in area.

Once the locations of the defects are obtained, the $+1/2$ defects are tracked
using the open source software Trackpy Allan  _et al._ (2021) using a
search_range value of 20 pixels. The trajectories thus obtained are further
filtered with a threshold of minimum three frames of survival.

Figure 4: Performance on simulation data as a function of window size.
$(1-R^{2})$ (blue circles) and average % error in coefficients (orange
triangles) for the (left) Q-tensor equation using the integral formulation and
(right) the flow equation using the weak formulation. The window sizes are
listed in pixels, while the simulation data used had a domain size of
$512\times 512\times 500$ pixels. Figure 5: Performance on simulation data as
a function of added noise. (left) $(1-R^{2})$ as a function of noise strength
for the orientation equation (blue circles), flow equation (orange diamonds)
and the flow equation in the weak form (green triangles). The first two are
obtained using the weak formulation. (right) The corresponding average % error
in the fitted coefficients. For the points absent on the plot, the obtained
optimal model did not match with the ground truth. A $5\times 5\times 5$
volume was used for the integral formulation, whereas a $205\times 205\times
65$ window was used for the weak formulation.

### Small field of view PolScope data

An additional dataset was taken using a combination of LC-PolScope microscopy
and dilute-labeled fluorescent MTs. With LC-PolScope microscopy, the
orientation field is obtained from birefringence information of polarized
light passing through the MT filaments that make up the nematic layer
Oldenbourg (2005) DeCamp  _et al._ (2015). This allows the orientation field
to be measured on MTs that are not fluorescent. However, a small fraction of
MTs in the sample were fluorescently labeled, and wide-field epifluorescence
images were acquired simultaneously with the birefringence data. We believe
that with wide-field microscopy, PIV on dilute-labeled MTs is more accurate
than on fully-labeled MTs, due to the difficulty of detecting velocity in the
direction of the elongated MT bundles.

This sample was prepared as detailed in Appendix C at 1.4 mM ATP with two
small differences: the chamber was created with double-sided tape and the
hydrophobic slide was made with Aquapel. Additionally, the proteins used in
this sample were from different preparations than the rest of the data in the
paper. The LC-PolScope sample was imaged on a Nikon Ti Eclipse with Andor Neo
camera, 20$\times$ 0.75 NA objective and 2 s frame interval.

Using the integral formulation with a window size of $5\times 5\times 5$
pixels, we obtain the equation

$\displaystyle\partial_{t}Q_{xx}=$ $\displaystyle-(1.07\
u_{x}\partial_{x}+1.08\ u_{y}\partial_{y})\ Q_{xx}$
$\displaystyle-(1.03\partial_{x}u_{y}-1.04\partial_{y}u_{x})\ Q_{xy}$
$\displaystyle+0.99\partial_{x}u_{x}$
$\displaystyle-2Q_{xx}\\{2.01Q_{xx}\partial_{x}u_{x}$
$\displaystyle\quad+Q_{xy}(1.00\partial_{x}u_{y}+1.05\partial_{y}u_{x})\\}$

with a high $R^{2}$ value of 0.97. This provides further strong evidence for
the discovered model. The flow analysis also yields a model similar to (5),
but with a low $R^{2}$, due to the limitation on the window size due to the
small field of view (see Fig. 4).

Figure 6: Identifying Q-tensor Equation using PolScope measurements.
Optimality curve for the $Q_{xx}$ equation for a small field-of-view dataset
in which the orientation is measured using PolScope microscopy. The $R^{2}$ at
the shoulder is $\sim 0.97$.

## References

  * Marchetti  _et al._ (2013) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, R. A. Simha, and M. Curie, Reviews of Modern Physics 85, 1143 (2013), arXiv:1207.2929 .
  * Ramaswamy (2010) S. Ramaswamy, Annual Review of Condensed Matter Physics 1, 323 (2010).
  * Toner  _et al._ (2005) J. Toner, Y. Tu, and S. Ramaswamy, Annals of Physics 318, 170 (2005).
  * Narayan  _et al._ (2007) V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 105 (2007).
  * Wensink  _et al._ (2012) H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Lowen, and J. M. Yeomans, Proceedings of the National Academy of Sciences 109, 14308 (2012).
  * Zhou  _et al._ (2014) S. Zhou, A. Sokolov, O. D. Lavrentovich, and I. S. Aranson, Proceedings of the National Academy of Sciences 111, 1265 (2014).
  * Duclos  _et al._ (2017) G. Duclos, C. Erlenkämper, J. F. Joanny, and P. Silberzan, Nature Physics 13, 58 (2017), iSBN: 1745-2473 1745-2481.
  * Kawaguchi  _et al._ (2017) K. Kawaguchi, R. Kageyama, and M. Sano, Nature 545, 327 (2017).
  * Kumar  _et al._ (2018) N. Kumar, R. Zhang, J. J. de Pablo, and M. L. Gardel, Science Advances 4, eaat7779 (2018).
  * Sanchez  _et al._ (2012) T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Nature 491, 431 (2012), arXiv:1301.1122 .
  * DeCamp  _et al._ (2015) S. J. DeCamp, G. S. Redner, A. Baskaran, M. F. Hagan, and Z. Dogic, Nat. Mater. 32, 1 (2015).
  * Doostmohammadi  _et al._ (2018) A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Nature Communications 9 (2018), 10.1038/s41467-018-05666-8.
  * Giomi  _et al._ (2013) L. Giomi, M. J. Bowick, X. Ma, and M. C. Marchetti, Physical Review Letters 110, 1 (2013).
  * Giomi  _et al._ (2014) L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek, and M. Cristina Marchetti, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20130365 (2014).
  * Doostmohammadi  _et al._ (2017) A. Doostmohammadi, T. N. Shendruk, K. Thijssen, and J. M. Yeomans, Nature Communications 8, 1 (2017).
  * Oza and Dunkel (2016) A. U. Oza and J. Dunkel, New Journal of Physics 18, 1 (2016).
  * Cortese  _et al._ (2018) D. Cortese, J. Eggers, and T. B. Liverpool, Physical Review E 97, 022704 (2018).
  * Shendruk  _et al._ (2018) T. N. Shendruk, K. Thijssen, J. M. Yeomans, and A. Doostmohammadi, Physical Review E 98, 010601 (2018), arXiv:1803.02093 .
  * Thampi  _et al._ (2014) S. P. Thampi, R. Golestanian, and J. M. Yeomans, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 90, 1 (2014).
  * Giomi (2015) L. Giomi, Physical Review X 5 (2015), 10.1103/PhysRevX.5.031003, arXiv:1409.1555 .
  * Lemma  _et al._ (2019) L. M. Lemma, S. J. DeCamp, Z. You, L. Giomi, and Z. Dogic, Soft Matter 15, 3264 (2019).
  * Shendruk  _et al._ (2017) T. N. Shendruk, A. Doostmohammadi, K. Thijssen, and J. M. Yeomans, Soft Matter 13, 3853 (2017), arXiv:1703.01531 .
  * Norton  _et al._ (2018) M. M. Norton, A. Baskaran, A. Opathalage, B. Langeslay, S. Fraden, A. Baskaran, and M. F. Hagan, Physical Review E 97, 1 (2018).
  * Gao  _et al._ (2017) T. Gao, M. D. Betterton, A.-S. Jhang, and M. J. Shelley, arXiv.org , 1 (2017), arXiv:1703.00969 .
  * Zhang  _et al._ (2016) R. Zhang, Y. Zhou, M. Rahimi, and J. J. de Pablo, Nature Communications 7, 13483 (2016).
  * Ellis  _et al._ (2017) P. W. Ellis, L. Giomi, D. J. G. Pearce, Y.-W. Chang, G. Goldsztein, and A. Fernandez-Nieves, Nature Physics 14, 85 (2017).
  * Alaimo  _et al._ (2017) F. Alaimo, C. Köhler, and A. Voigt, Scientific Reports 7, 1 (2017).
  * Brunton  _et al._ (2016) S. L. Brunton, J. L. Proctor, J. N. Kutz, and W. Bialek, Proceedings of the National Academy of Sciences of the United States of America 113, 3932 (2016), arXiv:1509.03580 .
  * Rudy  _et al._ (2017) S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Science Advances 3, 1 (2017), arXiv:1609.06401 .
  * Maddu  _et al._ (2022) S. Maddu, Q. Vagne, and I. F. Sbalzarini, arXiv:2201.08623 [cond-mat, physics:physics] (2022), arXiv:2201.08623 [cond-mat, physics:physics] .
  * Supekar  _et al._ (2021) R. Supekar, B. Song, A. Hastewell, A. Mietke, and J. Dunkel, arXiv:2101.06568 [cond-mat, physics:physics] (2021), arXiv:2101.06568 [cond-mat, physics:physics] .
  * Alves and Fiuza (2020) E. P. Alves and F. Fiuza, arXiv:2011.01927 [astro-ph, physics:physics] (2020), arXiv: 2011.01927.
  * Reinbold and Grigoriev (2019) P. A. Reinbold and R. O. Grigoriev, Physical Review E 100, 1 (2019), arXiv:1904.04314 .
  * Reinbold  _et al._ (2020) P. A. Reinbold, D. R. Gurevich, and R. O. Grigoriev, Physical Review E 101, 1 (2020), arXiv:1911.03365 .
  * Reinbold  _et al._ (2021) P. A. K. Reinbold, L. M. Kageorge, M. F. Schatz, and R. O. Grigoriev, Nature Communications 12, 3219 (2021).
  * Beris and Edwards (1994) A. N. Beris and B. J. Edwards,  _Thermodynamics of flowing systems with internal microstructure_ (Oxford University Press, 1994).
  * Lemma  _et al._ (2021) L. M. Lemma, M. M. Norton, A. M. Tayar, S. J. DeCamp, S. A. Aghvami, S. Fraden, M. F. Hagan, and Z. Dogic, Physical Review Letters 127, 148001 (2021).
  * Varghese  _et al._ (2020) M. Varghese, A. Baskaran, M. F. Hagan, and A. Baskaran, Physical Review Letters 125, 268003 (2020).
  * Zhou  _et al._ (2021) Z. Zhou, C. Joshi, R. Liu, M. M. Norton, L. Lemma, Z. Dogic, M. F. Hagan, S. Fraden, and P. Hong, Soft Matter 17, 738 (2021).
  * Duclos  _et al._ (2020) G. Duclos, R. Adkins, D. Banerjee, M. S. E. Peterson, M. Varghese, I. Kolvin, A. Baskaran, R. A. Pelcovits, T. R. Powers, A. Baskaran, F. Toschi, M. F. Hagan, S. J. Streichan, V. Vitelli, D. A. Beller, and Z. Dogic, Science 367, 1120 (2020).
  * Giomi  _et al._ (2012) L. Giomi, L. Mahadevan, B. Chakraborty, and M. F. Hagan, Nonlinearity 25, 2245 (2012), arXiv:arXiv:1110.4338v1 .
  * Giomi and Desimone (2014) L. Giomi and A. Desimone, Physical Review Letters 112, 1 (2014), _eprint: 1310.1908.
  * Chandragiri  _et al._ (2019) S. Chandragiri, A. Doostmohammadi, J. M. Yeomans, and S. P. Thampi, Soft Matter 15, 1597 (2019), arXiv:1901.06468 .
  * Chandragiri  _et al._ (2020) S. Chandragiri, A. Doostmohammadi, J. M. Yeomans, and S. P. Thampi, Physical Review Letters 125, 148002 (2020).
  * Thampi  _et al._ (2013) S. P. Thampi, R. Golestanian, and J. M. Yeomans, Physical Review Letters 111, 118101 (2013).
  * Thampi  _et al._ (2015) S. P. Thampi, A. Doostmohammadi, R. Golestanian, and J. M. Yeomans, EPL (Europhysics Letters) 112, 28004 (2015).
  * Thampi  _et al._ (2016) S. P. Thampi, A. Doostmohammadi, T. N. Shendruk, R. Golestanian, and J. M. Yeomans, Science Advances 2, e1501854 (2016).
  * Koch and Wilczek (2021) C.-M. Koch and M. Wilczek, Physical Review Letters 127, 268005 (2021).
  * Aditi Simha and Ramaswamy (2002) R. Aditi Simha and S. Ramaswamy, Physical Review Letters 89, 058101 (2002).
  * Maitra  _et al._ (2018) A. Maitra, P. Srivastava, M. Cristina Marchetti, J. S. Lintuvuori, S. Ramaswamy, and M. Lenz, Proceedings of the National Academy of Sciences of the United States of America 115, 6934 (2018).
  * Note (1) A higher accuracy result for the Q-tensor equation was obtained using a dataset with a small field of view with a few defects per frame (see Appendix Section E). That dataset however doesn’t give a good $R^{2}$ for the flow due to the limitation on the window size.
  * Note (2) We can force the framework to estimate dimensional parameters by constraining the regression procedure to include specific terms, while performing sparse regression on the remaining terms. For example, by forcing a term $K\nabla^{2}Q_{xx}$ (see Eq. (2)), we obtain a value for the elastic modulus of $K\sim$ $1\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}\text{\,}{\mathrm{s}}^{-1}$. However, because this term has a negligible contribution to the dynamics of $\mathbf{Q}$, the quantitative accuracy of this estimate may be limited.
  * Giomi  _et al._ (2011) L. Giomi, L. Mahadevan, B. Chakraborty, and M. F. Hagan, Physical Review Letters 106, 2 (2011).
  * Note (3) Independent measurements of the viscosity such as reported in Guillamat  _et al._ (2016) can be then used to estimate $\alpha$, the strength of the active force.
  * Note (4) For numerical stability, we include the $K\nabla^{2}\mathbf{Q}$ term in the $\mathbf{Q}$ equation with $K=1$.
  * Note (5) We re-scale the defect density so that the simulation value for $\alpha/\eta=0.3$ falls on the linear fit of the experimental values.
  * de and Prost (1993) G. P. G. de and J. Prost,  _The physics of liquid crystals_ (Clarendon Press, 1993).
  * Gao  _et al._ (2015) T. Gao, R. Blackwell, M. A. Glaser, M. D. Betterton, and M. J. Shelley, Physical Review Letters 114, 1 (2015), arXiv:1401.8059 .
  * Hemingway  _et al._ (2016) E. J. Hemingway, P. Mishra, M. C. Marchetti, and S. M. Fielding, Soft Matter 12, 7943 (2016).
  * Ngo  _et al._ (2014) S. Ngo, A. Peshkov, I. S. Aranson, E. Bertin, F. Ginelli, and H. Chaté, Physical Review Letters 113 (2014), 10.1103/PhysRevLett.113.038302.
  * Note (6) The second order term — $\lambda\\{\mathbf{Q}\cdot\mathbf{E}+\mathbf{E}\cdot\mathbf{Q}-2/d\Tr(\mathbf{Q}\cdot\mathbf{E})\mathbf{I}\\}$ Doostmohammadi  _et al._ (2018); Varghese  _et al._ (2020)—vanishes identically in 2D.
  * Zhao and Wang (2016) J. Zhao and Q. Wang, Journal of Scientific Computing 68, 1241 (2016).
  * Vanka (1986) S. P. Vanka, Journal of Computational Physics 65, 138 (1986).
  * Hermle (2017) R. Hermle, arXiv:1712.02030 [cs, math] (2017), arXiv:1712.02030 [cs, math] .
  * Castoldi and Popov (2003) M. Castoldi and A. V. Popov, Protein Expression and Purification 32, 83 (2003).
  * Hyman  _et al._ (1991) A. Hyman, D. Drechsel, D. Kellogg, S. Salser, K. Sawin, P. Steffen, L. Wordeman, and T. Mitchison, in  _Methods in Enzymology_, Molecular Motors and the Cytoskeleton, Vol. 196 (Academic Press, 1991) pp. 478–485.
  * Subramanian and Gelles (2007) R. Subramanian and J. Gelles, Journal of General Physiology 130, 445 (2007).
  * Lau  _et al._ (2009) A. W. C. Lau, A. Prasad, and Z. Dogic, EPL (Europhysics Letters) 87, 48006 (2009).
  * Rezakhaniha  _et al._ (2012) R. Rezakhaniha, A. Agianniotis, J. T. C. Schrauwen, A. Griffa, D. Sage, C. V. C. Bouten, F. N. van de Vosse, M. Unser, and N. Stergiopulos, Biomechanics and Modeling in Mechanobiology 11, 461 (2012).
  * Garcia (2011) D. Garcia, Experiments in Fluids 50, 1247 (2011).
  * Kamien (2002) R. D. Kamien, Reviews of Modern Physics 74, 953 (2002).
  * Allan  _et al._ (2021) D. B. Allan, T. Caswell, N. C. Keim, C. M. van der Wel, and R. W. Verweij, “Soft-matter/trackpy: Trackpy v0.5.0,” Zenodo (2021).
  * Oldenbourg (2005) R. Oldenbourg, Live cell imaging: a laboratory manual , 205 (2005).
  * Guillamat  _et al._ (2016) P. Guillamat, J. Ignés-Mullol, S. Shankar, M. C. Marchetti, and F. Sagués, Physical Review E 94, 1 (2016), arXiv:1606.05764 .