1]Center for Soft Matter Research, Department of Physics, New York University,
New York, NY,10003, USA

# Model chromatin flows: numerical analysis of linear and nonlinear
hydrodynamics inside a sphere

Iraj Eshghi<EMAIL_ADDRESS>Alexandra Zidovska<EMAIL_ADDRESS>Alexander Y.
Grosberg<EMAIL_ADDRESS>[

###### Abstract

We solve a hydrodynamic model of active chromatin dynamics, within a confined
geometry simulating the cell nucleus. Using both analytical and numerical
methods, we describe the behavior of the chromatin polymer driven by the
activity of motors having polar symmetry, both in the linear response regime
as well as in the long-term, fully nonlinear regime of the flows. The
introduction of a boundary induces a particular geometry in the flows of
chromatin, which we describe using vector spherical harmonics, a tool which
greatly simplifies both our analytical and numerical approaches. We find that
the long-term behavior of this model in confinement is dominated by steady,
transverse flows of chromatin which circulate around the spherical domain.
These circulating flows are found to be robust to perturbations, and their
characteristic size is set by the size of the domain. This gives us further
insight into active chromatin dynamics in the cell nucleus, and provides a
foundation for development of further, more complex models of active chromatin
dynamics.

##  1 Introduction

The cell nucleus houses the genome, containing genetic information needed for
the cell’s life alberts2017molecular . This information is encoded into a long
DNA molecule, which in cells forms a complex with histone proteins, the
chromatin fiber van2012chromatin . The cell nucleus encompasses chromatin as
well as a variety of molecules such as proteins and RNA alberts2017molecular ;
spector1993macromolecular , and is maintained out-of-equilibrium by a large
number of active processes e.g. transcription, replication, chromatin
remodeling and DNA repair rippe2007dynamic ; zidovska2020rich . Active
processes affect both the dynamics hubner2010chromatin ; fierz2019biophysics ;
zidovska2020self  and organization of the genome hubner2013chromatin ;
bickmore2013genome ; mirny2022mechanisms . Examples include changes in:
chromatin mobility with transcription soutoglou2007mobility ;
lanctot2007dynamic , local chromatin packing and dynamics upon DNA damage
mine2013dna ; lebeaupin2015chromatin ; eaton2020structural , occasional
directed motion of chromosomal loci marshall1997interphase ; levi2005chromatin
; chuang2006long , or the active extrusion of chromatin loops by cohesin
mirny2021keeping ; davidson2021genome . In particular, active processes were
shown to lead to the formation of micron-scale chromatin domains, within which
chromatin moves coherently zidovska2013micron .

The first theoretical model developed to explore the origins of the chromatin
coherent motions was based on two-fluid hydrodynamics bruinsma2014chromatin ,
it was followed by explicit hydrodynamic computational models
saintillan2018extensile ; Mahajan2022hetero . These papers led to the
important insight that local active forces mediated through long-range
hydrodynamic interactions can lead to large-scale motions of chromatin.
Separately, several computational models were developed which do not include
hydrodynamic interactions, but reproduce the phenomenology of coherent motions
at the expense of introducing artificial long-range forces between sections of
the polymer shi2018interphase ; liu2018chain ; di2018anomalous .

An analytical description of chromatin hydrodynamics was developed in the
works bruinsma2014chromatin ; eshghi2022symmetry ; eshghi2023alignment  based
on the two-fluid model doi1992dynamic  and the ideas of statistical physics of
active systems marchetti2013hydrodynamics ; vicsek1995novel ; toner1995long ;
AdarJoanny2021 . The most significant flows were obtained under the assumption
that active motors driving are such that each of them exerts one force on the
polymer and equal but opposite force on the solvent. Indeed typical active
motors operating on the genome can be viewed as such, for instance, RNA
polymerase II which is a common active enzyme in the nucleus
kimura1999quantitation . Therefore, in both our previous work
eshghi2023alignment  as well as in the present paper, we consider our theory
in the most general phenomenological form, not specifying the nature of motors
beyond the fact that they have a polar symmetry. It is because of polar
symmetry that they act on the relative velocity of the solvent past polymer
like force monopoles and not dipoles, thus generating very significant flows.
If the number and activity of motors exceeds a critical threshold, we found
that they spontaneously form an ordered polar phase which actively pumps
chromatin and solvent through one another. We successfully described this
spontaneous ordering through a polar order parameter, and identified the value
of this critical threshold as a function of model parameters along with the
critical exponents near the transition eshghi2023alignment .

These results, however, were found with the assumption of an infinite
boundless medium. While this made the analysis more straightforward, a more
accurate description of the situation experienced by the chromatin polymer
would account for the finite size of the cell nucleus. In this work, we will
study the dynamics of our model in a confined spherical geometry, near its
critical point. We show that the critical value of force and density of motors
is shifted by an amount dependent on model parameters and system size.
Furthermore, we find that the length scale of the modes excited near the
transition point is set by the confinement size, consistently with
expectations from the equilibrium theory of second-order phase transitions.

##  2 The model, confining geometry, and equations of motion

Our goal in the present work is to analyze previously derived equations of
hydrodynamic motion of chromatin eshghi2023alignment , both linearized and
nonlinear, when confined in a spherical domain. To make this work self-
contained, and to set up the notations, we start by re-stating the primary
equations of motion of this model.

###  2.1 Equations of motion

We consider two mutually permeating fluids: a polymer and a solvent. Their
velocities are denoted by the fields
$\mathbf{v}^{\mathrm{p}},\;\mathbf{v}^{\mathrm{s}}$ respectively, although in
our calculations we will use the following two linear combinations of the
velocity fields:
$\mathbf{w}=\mathbf{v}^{\mathrm{p}}-\mathbf{v}^{\mathrm{s}},\;\mathbf{u}=(\eta^{\mathrm{p}}\mathbf{v}^{\mathrm{p}}+\eta^{\mathrm{s}}\mathbf{v}^{\mathrm{s}})/(\eta^{\mathrm{p}}+\eta^{\mathrm{s}})$.
The two fluids, flowing past another, experience a friction per unit volume
$\zeta$. The volume fraction of the polymer is $\phi(\mathbf{r})$, and the
two-fluid combination is assumed to be incompressible, so the volume fraction
of solvent is $1-\phi(\mathbf{r})$. Both fluids are assumed to experience
viscous dissipation upon shear, with respective viscosities
$\eta^{\mathrm{p}},\;\eta^{\mathrm{s}}$. The viscosity of the chromatin
polymer is known to have a frequency dependence eshghi2021interphase , but
here we will assume that it is simply Newtonian. The polymer experiences an
osmotic pressure $\Pi(\phi)$, which is only a function of $\phi$ as it is
assumed to equilibrate quickly and locally. We consider a regime where the
polymer density deviates weakly from its mean value $\phi_{0}$, allowing us to
linearize the equations of motion around that point:
$\phi=\phi_{0}+\delta\phi$. As a result the osmotic pressure can be written as
$\Pi=\Pi_{0}+K\delta\phi$, where $K$ is the osmotic modulus of the polymer.
This gives us the following equations for the fluids

$\displaystyle\begin{split}&\zeta\left(\frac{1}{\eta^{\mathrm{p}}}+\frac{1}{\eta^{\mathrm{s}}}\right)\mathbf{w}=\nabla^{2}\mathbf{w}+\left(\frac{\mathbf{F}^{\mathrm{p}}}{\eta^{\mathrm{p}}}-\frac{\mathbf{F}^{\mathrm{s}}}{\eta^{\mathrm{s}}}\right)\\\
&\ \ \ \ \ \ \ \ \ \ \
+\left(\frac{1-\phi_{0}}{\eta^{\mathrm{s}}}-\frac{\phi_{0}}{\eta^{\mathrm{p}}}\right)\nabla
P-\frac{K}{\eta^{\mathrm{p}}}\nabla\delta\phi\ ,\end{split}$ (1a)
$\displaystyle\nabla\cdot\mathbf{w}=-\partial_{t}\delta\phi\left(\frac{1}{\phi_{0}}+\frac{1}{1-\phi_{0}}\right)\
,$ (1b)
$\displaystyle\left(\eta^{\mathrm{p}}+\eta^{\mathrm{s}}\right)\nabla^{2}\mathbf{u}=K\nabla\delta\phi+\nabla
P+\mathbf{F}^{\mathrm{p}}+\mathbf{F}^{\mathrm{s}}\ ,$ (1c)
$\displaystyle\left(\eta^{\mathrm{p}}+\eta^{\mathrm{s}}\right)\nabla\cdot\mathbf{u}=-\partial_{t}\delta\phi\left(\frac{\eta^{\mathrm{p}}}{\phi_{0}}-\frac{\eta^{\mathrm{s}}}{1-\phi_{0}}\right)\
,$ (1d)

where $\mathbf{F}^{\mathrm{p}},\;\mathbf{F}^{\mathrm{s}}$ are the forces
densities generated by the motors, and $P$ is the hydrostatic pressure induced
by total incompressibility.

As stated in the introduction, we assume that every active motor exerts a
force on the polymer and equal and opposite force on the solvent. In the
coarse-grained description, assuming there are some $\rho$ active motors per
unit volume, we can write

$\mathbf{F}^{p}(r)=-\mathbf{F}^{s}(r)=\rho f\mathbf{m}(r)\ ,$ (2)

where $f$ is the force produced by a single motor, while $\mathbf{m}(r)$ is
the average vector of orientation of motors located around point $\mathbf{r}$.
According to equations (1a) and (1c), field $\mathbf{w}$ is driven by force
monopole density $f\rho\mathbf{m}$, while field $\mathbf{u}$ in this
approximation is not driven at all.

Vector field $\mathbf{m}(\mathbf{r})$ plays an important role in our theory,
and it is worth several comments. First, the very possibility to define vector
$\mathbf{m}$ is because motors we consider have polar symmetry, they have two
different sides, one attached to the polymer and another exposed to the
solvent. Second, vector $\mathbf{m}$, as the average of the unit orientation
vectors of individual motors, has magnitude that is always smaller than unity:
$|\mathbf{m}|\leq 1$. Third, this vector naturally serves as an order
parameter of polarization ordering phase transition predicted and described by
our theory.

###  2.2 Linear Dynamics

Given the polar symmetry of motors, each of them is subject to a torque
whenever there is a flow of the solvent relative to the polymer, i.e., when
$\mathbf{w}\neq 0$. If drive is weak (because motors are either weak or not
numerous enough), then this torque is also weak and the distribution of motor
orientations is nearly isotropic. In this case, the dynamics of average
orientation is described by

$\tau\partial_{t}\mathbf{m}=-2\mathbf{m}+\frac{2\tau}{3a}\mathbf{w}\;,$ (3)

where we have introduced the reorientation time of the dipoles
$\tau=\gamma/T$. $\gamma$ is the rotational friction coefficient of the motor,
$T$ is the ambient (potentially effective) temperature, and $a$ is the size of
the motors, assumed to be comparable to (or smaller than) the mesh size of the
polymer. We can now combine equations (3) and (1) into one equation for
$\mathbf{w}$ and $\mathbf{m}$.

$\begin{split}&\left(1+\lambda^{2}\nabla\times\nabla\times-\lambda_{s}^{2}\nabla\nabla\cdot\right)\tau\partial_{t}\mathbf{w}\\\
&-2\lambda_{d}^{2}\nabla\nabla\cdot\mathbf{w}_{\omega}=\frac{f\rho\tau}{\zeta}\partial_{t}\mathbf{m}_{\omega}\;,\end{split}$
(4)

where we have defined three length scales:

$\lambda^{2}=\frac{\eta^{\mathrm{s}}}{\zeta},\quad\lambda_{s}^{2}\simeq\frac{\eta^{\mathrm{p}}\left(1-\phi_{0}\right)^{2}}{\zeta},\
\ \text{and}\ \
\lambda_{d}^{2}\simeq\frac{K\phi_{0}\left(1-\phi_{0}\right)^{2}\tau}{2\zeta}\
.$ (5)

The first of these, $\lambda$, is naturally identified as the mesh size of the
polymer. Second, significantly larger length scale is $\lambda_{s}$, the
screening length of hydrodynamic interactions in the two-fluid system.
Finally, the third length scale, $\lambda_{d}$, characterizes the interplay of
the two fluid system with motors, namely, it is the typical distance of
cooperative diffusion by the polymer during motor reorientation time $\tau$;
of course, this motion can be thought of as driven by osmotic elasticity of
the polymer (described by $K$) and opposed by friction (described by $\zeta$).
Since we are treating chromatin in a continuum approximation, the model we are
working with is only relevant on length scales larger than the mesh size
$\lambda$.

Since all of the equations of motion in this regime are linear, we can
simultaneously decompose the fields $\mathbf{u},\;\mathbf{w},\;\mathbf{m}$
into their divergence-free (transverse, $\perp$) and curl-free (longitudinal,
$\parallel$) components. Separating these into their respective equations of
motion, we obtain

$\displaystyle\qquad\quad\left(1-\lambda_{s}^{2}\nabla^{2}\right)\mathbf{w}_{\perp}$
$\displaystyle=\frac{f\rho}{\zeta}\mathbf{m}_{\perp}$ (6a)
$\displaystyle\left(1-\lambda_{s}^{2}\nabla^{2}\right)\tau\partial_{t}\mathbf{w}_{\parallel}$
$\displaystyle=2\lambda_{d}^{2}\nabla^{2}\mathbf{w}_{\parallel}+\frac{f\rho\tau}{\zeta}\partial_{t}\mathbf{m}_{\parallel}\;.$
(6b)

Our previous study allowed us to identify the value of the force $f$ above
which an ordered phase for $\mathbf{m}$ spontaneously develops. In an infinite
domain, this occurs at $f=3a\zeta/\rho\tau$, where the velocity generated by
the force dipoles $f\rho/\zeta$ equals the velocity needed to orient them
$a/\tau$. This allows us to identify the critical parameter

$\epsilon=\frac{f\rho\tau}{3a\zeta}-1\;.$ (7)

When $\epsilon>0$, the zero modes of the velocity and orientation fields
become unstable. This all changes, however, with the addition of boundaries,
which set a maximal length scale for any dynamics for the two fluids. Our goal
in this paper is to refine the conclusions of our previous work in the finite
domain context.

###  2.3 Nonlinear Dynamics

Once an instability develops, the polar order parameter grows in magnitude
exponentially. Eventually, nonlinear effects inevitably kick in. In
particular, they guarantee that $|\mathbf{m}|\leq 1$. The precise evolution of
$\mathbf{m}$ in the nonlinear regime is complex, as it couples to all the
other moments of the distribution of motor orientations. However, in our
previous work eshghi2023alignment  we derived an approximate equation of
motion for $\mathbf{m}$ including only a second-order nonlinearity in
$\mathbf{w}$:

$\begin{split}\tau\partial_{t}\mathbf{m}&=2(\mathbf{m}_{eq}(\mathbf{w})-\mathbf{m})\\\
\mathbf{m}_{eq}&=\frac{\mathbf{w}\tau}{3a}\left(1-\frac{3}{5}\left(\frac{\mathbf{w}\tau}{3a}\right)^{2}\right)\;.\end{split}$
(8)

This is just one of a long list of possible nonlinearities that may be
considered. We choose to consider it above all others as it naturally controls
the amplitude of the unstable modes, and without the addition of any further
nonlinearities leads to stable long-term dynamics for this system. While
$\mathbf{m}$ obeys equation (8), the velocity field $\mathbf{w}$ still follows
equation (4). This system cannot be solved analytically in general, so we will
turn to a numerical method in the following section.

###  2.4 Boundary Conditions

Even though the nucleus usually looks like an ellipsoid, we reduce the
complexity of the problem by modeling it as a sphere of radius $R$.

We assume no-slip boundary conditions for both velocity fields
$\mathbf{v}^{\mathrm{p}},\;\mathbf{v}^{\mathrm{s}}$. In the case of the
polymer, this is justified by the tethering of chromatin to the boudary by
LINC complexes crisp2006coupling , while no-slip boundary conditions for the
solvent are standard for viscous fluids batchelor1967introduction

$\mathbf{v}^{\mathrm{p}}_{\mathrm{tangential}}\big{|}_{r=R}=\mathbf{v}^{\mathrm{s}}_{\mathrm{tangential}}\big{|}_{r=R}=0\;.$
(9)

Turning now to normal components of velocity, we assume that no permeation of
the boundary is possible by either polymer or solvent. We make this assumption
despite the fact that the nuclear membrane is porous and lets some small
molecules through alberts2017molecular . We do so not only for simplicity, but
also because over the seconds-timescale that we are interested in the volume
of the nucleus is conserved chu2017origin , thus the net flux of material
through it is small. Finally, we treat the nuclear boundary as rigid and
ignore its fluctuations, which are another source of potentially interesting
effects chu2017origin . All these simplifications result in the following
boundary condition:

$\left[\phi\mathbf{v}^{\mathrm{p}}+(1-\phi)\mathbf{v}^{\mathrm{s}}\right]_{\mathrm{normal}}|_{r=R}=0\;.$
(10)

The meaning of this condition is simple: although neither component can go
through the membrane, any one component, either polymer or solvent, can be
approaching the membrane with some non-zero velocity provided that the other
component at the same time departs from the wall, such that the exchange
between them does not change volume.

##  3 Results

The spatial behavior of the linearized model (6a,6b) is entirely contained in
Laplacian operators. We will use the vector eigenfunction spectrum of the
Laplacian to construct exact solutions for this model in the linear response
regime. Given the restriction of the solutions to a closed bounded domain, the
spectrum is discrete.

Once we turn to the nonlinear regime, it becomes impossible to solve the
equations of motion exactly using these basis functions. However, we have
constructed a numerical method for the full nonlinear model which exploits the
basis we develop for the linearized version. This makes three-dimensional
solutions of the full nonlinear PDE possible with minimal computational
complexity.

###  3.1 Linearized model: an analytical solution

Figure 1: First excited longitudinal and transverse modes, shown along a
vertical slice of the spherical domain. Since the vector dependence of these
first modes is simple, we choose to plot only one component of the resulting
vector fields, since the other components are simply $0$. In the case of the
first longitudinal mode the flow is spherically symmetric, while the first
transverse mode is axially symmetric about the $\hat{\mathbf{z}}$ axis.

We begin by finding solutions in the linear response regime, where the
distribution of motor orientations deviates weakly from uniform. This allows
us to expand the flow fields in a basis which simplifies the equation of
motion. If we consider curl-free ($\mathbf{w}_{\parallel}$) and divergence-
free ($\mathbf{w}_{\perp}$) components of the velocity separately, then the
dynamics of these modes separate and simplify, as shown in equations (6).

The natural choice of basis functions to solve this system of equations needs
to have two properties: they need to be eigenfunctions of the Laplace operator
in spherical coordinates, to simplify the spatial dependence of the equations
of motion, and they must be vectors to preserve the symmetry of the velocity
fields. Such functions are already known, and are referred to as vector
spherical harmonics (VSH) mcc1981morse . They are constructed based on the
well-known scalar spherical harmonics $Y_{lm}(\theta,\phi)$. There are
several, closely related, possible definitions for VSH. We will use the
following in this paper:

$\begin{split}\mathbf{Y}_{lm}&=Y_{lm}\hat{\mathbf{r}},\\\
\mathbf{\Psi}_{lm}&=r\nabla Y_{lm},\\\
\mathbf{\Phi}_{lm}&=\mathbf{r}\times\nabla Y_{lm}\;.\end{split}$ (11)

From these, we construct curl-free and divergence-free components of
$\mathbf{w}$. First, notice that for any scalar function $f(\mathbf{r})$, we
have $\nabla\cdot\left(f\mathbf{\Phi}_{lm}\right)=0$. Then recall that for a
vector-valued function to be curl-free, it suffices that it be the gradient of
a scalar. At every radial shell $r$, we expand the scalar function
$f(\mathbf{r})$ in spherical harmonics, with coefficients $f_{lm}(r)$. This
gives us

$\begin{split}\nabla
f(\mathbf{r})&=\sum_{lm}\nabla\left(f_{lm}(r)Y_{lm}\right)\\\
&=\sum_{lm}\left(\frac{\partial f_{lm}(r)}{\partial
r}\mathbf{Y}_{lm}+\frac{f_{lm}(r)}{r}\mathbf{\Psi}_{lm}\right)\;.\end{split}$
(12)

Therefore, we will solve equations (6) using the following expansions:

$\begin{split}\mathbf{w}_{\perp}&=\sum_{lm}a_{lm}(r,t)\mathbf{\Phi}_{lm};\\\
\mathbf{w}_{\parallel}&=\sum_{lm}\nabla\left(b_{lm}(r,t)Y_{lm}\right)\\\
&=\sum_{lm}\left(\frac{\partial b_{lm}(r,t)}{\partial
r}\mathbf{Y}_{lm}+\frac{b_{lm}(r,t)}{r}\mathbf{\Psi}_{lm}\right)\;.\end{split}$
(13)

To identify what spatial dependence $a_{lm},b_{lm}$ must have, we search for
eigenfunctions of the Laplacian which meet our boundary conditions. Here it is
important to notice that by construction, the longitudinal velocity field
$\mathbf{w}_{\parallel}$ automatically obeys condition (10) due to the
continuity of $\delta\phi$. Therefore, we only need to impose (9) in the case
of the longitudinal flows. As a result, we are searching for
$a_{lm}(r,t),b_{lm}(r,t)$ which satisfy

$\begin{split}&a_{lm}|_{r=R}=b_{lm}|_{r=R}=0\\\
&\nabla^{2}\left(a_{lm}\mathbf{\Phi}_{lm}\right)=-k^{2}a_{lm}\mathbf{\Phi}_{lm}\;,\\\
&\nabla^{2}\left(b_{lm}Y_{lm}\right)=-k^{2}b_{lm}Y_{lm}\;,\end{split}$ (14)

for some real number $k^{2}$. One such set of functions are spherical Bessel
functions $j_{l}(x)$, which have the property
$\nabla^{2}\left(j_{l}(kr)Y_{lm}\right)=-k^{2}j_{l}(kr)Y_{lm}$ in spherical
coordinates. We choose $k$ such that the modes meet the boundary conditions at
$r=R$. This results in the following expansion

$\begin{split}\mathbf{w}_{\perp}&=\sum_{lmn}a_{lmn}(t)j_{l}\left(\frac{\alpha_{ln}r}{R}\right)\mathbf{\Phi}_{lm}\;,\\\
\mathbf{w}_{\parallel}&=\sum_{lmn}{b}_{lmn}(t)\nabla\left(j_{l}\left(\frac{\alpha_{ln}r}{R}\right)Y_{lm}\right)\;,\\\
\end{split}$ (15)

where we have defined $\alpha_{ln}$ to be the $n$th zero of the $l$th order
spherical Bessel function. Equipped with the basis (15), it is now
straightforward to insert the expansion into the equations of motion (6), and
solve for the time dependence of the coefficients $a_{lmn},b_{lmn}$:

$\displaystyle\tau\dot{a}_{lmn}(t)=\frac{2\left(\epsilon-\lambda^{2}\left(\alpha_{ln}/R\right)^{2}\right)}{1+\lambda^{2}\left(\alpha_{ln}/R\right)^{2}}a_{lmn}(t)\;,$
(16a)
$\displaystyle\begin{split}&\tau^{2}\ddot{b}_{lmn}(t)-2\frac{\left(\epsilon-(\lambda_{d}^{2}+\lambda_{s}^{2})\left(\alpha_{ln}/R\right)^{2}\right)}{1+\lambda_{s}^{2}\left(\alpha_{ln}/R\right)^{2}}\tau\dot{b}_{lmn}(t)\\\
&\quad+4\frac{\lambda_{d}^{2}\left(\alpha_{ln}/R\right)^{2}}{1+\lambda_{s}^{2}\left(\alpha_{ln}/R\right)^{2}}b_{lmn}(t)=0\;.\end{split}$
(16b)

The time dependent dynamics of these coefficients and the corresponding modes
(15) is very similar to what we have found previously for an unbounded domain,
except of course the instability thresholds are significantly affected due the
effect of boundaries. In the case of the transverse modes, the solutions of
equation (16a) are simple exponentials in time, with the sign of the
coefficient on the right-hand-side of equation (16a) determining stability.
While at small $\epsilon$ the mode is stable, as soon as
$\epsilon>\lambda^{2}\alpha_{ln}^{2}/R^{2}$ the mode becomes unstable. As
expected, the finite-size effect is captured by the unitless ratio
$\lambda/R$. If the domain is a lot larger than the mesh size, which
corresponds to the thickness of the boundary layer needed to meet the no-slip
condition, then the dynamics are similar to that of an infinite domain.
However, if the domain is small, additional forcing is needed to excite these
modes, since the viscous drag at the boundary will significantly dampen their
motion.

The time-dependence of the longitudinal modes, as in the infinite domain case
(see eshghi2023alignment ), is that of a harmonic oscillator. The term in
equation (16b) which controls stability is the friction term, and it flips
sign when $\epsilon$ is sufficiently large resulting in an apparent negative
friction coefficient. As a result, the amplitude of the oscillations grow
exponentially in time. The necessary forcing to drive this instability is far
larger than the transverse case. Indeed, the critical value of $\epsilon$ in
this case is $(\lambda_{s}^{2}+\lambda_{d}^{2})\alpha_{ln}^{2}/R^{2}$, and
recall that $\lambda_{s}\gg\lambda$. Each of those terms represents a source
of dissipation which must be overcome. $\lambda_{s}$ reflects the energy
dissipation due to friction between polymer and solvent in these longitudinal
modes. This term is larger than in the case of the transverse modes (in that
case it was proportional to $\lambda$), because the transverse velocity of
polymer is far smaller than that of solvent due to force balance and the fact
that $\eta^{\mathrm{p}}\gg\eta^{\mathrm{s}}$. In the case of longitudinal
flows however, the polymer can compress and move faster, leading to increased
friction. Finally, $\lambda_{d}$ reflects the dissipation due to density
change of the polymer.

It is worth considering what the first excited modes are in this geometry. The
lowest-$\epsilon$ (easiest to excite) transverse mode is proportional to
$\mathbf{\Phi}_{10}=\sin(\theta)\hat{\mathbf{\phi}}$. In this state, the
chromatin is swirling around the $\hat{z}$ axis, with highest speed around the
equator and decaying to $0$ at the poles. In contrast, the first longitudinal
mode to be excited is proportional to
$\mathbf{Y}_{00}=\frac{1}{4\pi}\hat{\mathbf{r}}$. This describes a radially-
symmetric ”breathing” motion of the chromatin, moving in and out of the center
of the nucleus in an oscillatory manner. We show a slice of these first
excited modes in Figure 1.

In summary, restoring physical parameters, the condition for the active force
density is

$f\rho>\frac{3a\zeta T}{\gamma}+\frac{a\eta_{s}T\alpha_{11}^{2}}{\gamma
R^{2}}$ (17)

for the transverse modes (the lowest allowable mode in this case is $l=1$,
since $\mathbf{\Phi}_{00}=0$). In the above, we have $\alpha_{11}\simeq 4.5$.
For the longitudinal modes, the condition is

$f\rho>\frac{3a\zeta
T}{\gamma}+\frac{(1-\phi_{0})^{2}a\pi^{2}}{R^{2}}\left[\frac{3T\eta^{\mathrm{p}}}{\gamma}+\frac{3}{2}K\phi_{0}\right]\;,$
(18)

where we have inserted the value of $\alpha_{01}=\pi$. These are the same
expressions that we previously guessed eshghi2023alignment , except the
factors of $\alpha_{ln}^{2}$ were absent. We find that the addition of
boundaries induces necessary gradients in the velocity field which the sources
of activity must compete against. Furthermore, the necessary excess force is
far larger if the motors are to compress and pump the polymer in a finite
domain. Consistently with our previous results, these effects become less
significant the larger the domain size is compared to the length scales
$\lambda,\lambda_{s},\lambda_{d}$.

###  3.2 Nonlinear regime

Figure 2: Transverse flows are stable in the nonlinear regime. A: Amplitude of
the coefficients for the lowest transverse and longitudinal modes ($a_{10}$
and $b_{01}$ respectively), as a function of time when initialized in a pure
state with $a_{10}(0)=0.1$. Nonlinearities mix multiple transverse modes
together, but no longitudinal modes are excited. B: Amplitude of the same
modes as in A, but with the initial conditions set to random numbers of
magnitude $10^{-8}$, with $a_{10}(0)=0.1$. The transverse steady-state is
robust to this small perturbation.

Now that we have identified the linear stability of this system with
$\epsilon$ slightly above the critical point, we must consider the long-time
behavior once the instabilities have grown and saturated the order parameter
field $\mathbf{m}$. The nonlinear dynamics of this system, described by
equations (4,8), have to be investigated using a numerical method. We
developed such a method, which we describe in detail in the Appendix, section
B. To integrate the equations of motion, we convert the fields from a spatial
representation to the basis (15) at each time-step, then evolve the
coefficients, before transforming back to spatial representation to evaluate
the nonlinear time evolution for $\mathbf{m}$ shown above. We measured all
length scales relative to the confinement scale $R$ and chose the numerical
values $\lambda_{d}/R=10^{-1}$, $\lambda_{s}/R=10^{-2}$, and
$\lambda/R=10^{-3}$. This choice is consistent with our estimates that
$\lambda_{s}\gg\lambda$ and also with numerical values of average mesh size in
chromatin $\lambda\approx 70\,\mathrm{nm}$ and $R\approx 10\,\mathrm{\mu m}$.
We set $\epsilon=0.3$, above both thresholds of instability (17,18).

The first question we set to explore is whether the modes identified in the
linear response regime become stable once the nonlinearity becomes
significant. We first test the transverse modes, by starting the system in a
pure state $\mathbf{w}(t=0)=w_{0}j_{1}(r\alpha_{11}/R)\mathbf{\Phi}_{10}$,
with $w_{0}=0.1$ and allowing it to evolve for $100$ time-steps of step size
$\Delta t=0.5$. Eventually the dynamics settle into a steady-state consisting
of a superposition of transverse states, but no significant longitudinal
excitations are observed, as shown in Figure 2A. This result is robust to
small perturbations, which we test by initializing all other modes with
uniformly distributed random values ranging between $\pm 10^{-8}$, as shown in
Figure 2B.

We then perform the same analysis but with the longitudinal modes. We
initialize the system in a pure state
$\mathbf{w}(t=0)=w_{0}\nabla\left(j_{0}(r\alpha_{01}/R)Y_{00}\right)$ with
$w_{0}=0.1$. To better resolve the oscillatory dynamics in this case we lower
the time step to $\Delta t=0.25\tau$, and we evolve for $200$ steps. The
oscillations initially grow in magnitude until they settle at a maximum, at
which point their amplitude remains steady for all times we integrated, as
shown in Figure 3A. However, when we initialize the configuration with uniform
random numbers ranging between $\pm 10^{-8}$ for all other modes, after a
certain number of cycles the longitudinal oscillator gets taken over by the
transverse steady-state modes, as shown in Figure 3B. To see where the
evolution stabilizes we keep the system going for $2000$ time steps, and find
that while the root-mean-squared magnitude of transverse modes reaches
something close to a steady-state, individual modes such as the first
coefficient $a_{10}(t)$ still evolve over those timescales, shown in Figure
3B. To verify that they eventually do reach a steady-state, we increase the
time step to $\Delta t=0.5\tau$ and evolve for $10^{4}$ steps, finding that
the modes do eventually fully settle, albeit in a different steady-state each
time due to the random initial conditions. One example evolution is shown in
Figure 3C.

Figure 3: Longitudinal flows are unstable to perturbations in the nonlinear
regime. A: Amplitude of the coefficients for the lowest transverse and
longitudinal modes ($a_{10}$ and $b_{01}$ respectively), as a function of time
when initialized in a pure state with $b_{10}(0)=0.1$. Nonlinearities mix
longitudinal modes together, but the system remains free of transverse modes.
B: Amplitude of the same modes as in A, but with the initial conditions set to
random numbers of magnitude $10^{-8}$, with $b_{10}(0)=0.1$. In this case, the
transverse flows grow over time and eventually dominate the evolution of the
system, suppressing all oscillations. C: Same as in B, but evolved over even
longer timescales, to verify that the system does indeed reach a steady-state
eventually.

We test this stability further by allowing the longitudinal oscillations to
fully develop before introducing transverse flows. We initialize the system in
a pure longitudinal state as before, and allow it to evolve for $150$ time
steps with step size $\Delta t=0.25\tau$, which is equivalent to about three
full cycles of the oscillator. Then, we introduce a transverse kick by setting
the amplitude of the lowest transverse mode to a large value, $a_{10}=3$. This
briefly perturbs the oscillations as can be seen in the amplitude of $b_{01}$
right afterwards in Figure 4, but the system quickly settles back into its
previously established resonance. Evidently, once the longitudinal oscillatory
flows have settled in place, they are more robust to perturbations.

##  4 Discussion & Conclusions

In this paper we have analyzed solutions of the active two-fluid equations of
motion we derive in eshghi2023alignment , but in a confined domain. As this
model is intended to describe the active dynamics of the chromatin polymer and
its solvent nucleoplasm, the addition of the confinement allows us to study
the effect of the nuclear boundary on these dynamics. Specifically, we are
interested in the kinds of active flows which may take place in such a
restricted environment, and the ways in which the hierarchies of length scales
(mesh size, screening length, osmotic relaxation given a characteristic time)
may affect the stability of the dynamics. In our previous work
eshghi2023alignment , we performed estimates of the relevant length, time, and
velocity scales using experimental values and found the theory to be
consistent with experiments.

As estimated in that work, we find that the confinement length scale $R$ leads
to factors of $1/R^{2}$ everywhere where the Laplace operator is present,
resulting in confinement-dependent instability thresholds for both
longitudinal and transverse modes. In the limit $R\rightarrow\infty$, we
recover the unbounded results derived in eshghi2023alignment . Furthermore, we
find that the confinement size determines the size of the eigenmodes which get
excited, which will be an interesting result to compare to experimental data
on active chromatin flows.

Figure 4: Longitudinal oscillations are stable to late perturbations. The
system is prepared in the same manner as before, in a pure longitudinal state,
but at a prescribed time the transverse mode $a_{10}$ is arbitrarily increased
in magnitude to introduce a kick to teh system. For a short period the lowest
longitudinal mode amplitude (orange) is affected, but it quickly settles back
into its previous oscillatory behavior.

We go beyond linear stability analysis by numerically integrating the
equations of motion in a full three-dimensional domain. We find that the
decomposition into transverse and longitudinal modes, which we performed for
the sake of linear stability analysis, remains relevant in the nonlinear case.
In particular we find that if initialized in a pure transverse or longitudinal
state, the system will remain in that class of states (although the
nonlinearity will mix transverse modes among one another and equivalently for
longitudinal modes). Transverse modes are found to be stable under small
perturbations, while longitudinal modes are stable under some but not others.
It will be interesting, in a later study, to analyze in detail the basins of
attraction and regimes of stability for all of these types of flows in the
nonlinear regime. In sum, we find that the minimal active two-fluid model we
have developed, exhibits rich and nontrivial behavior.

It is worth asking whether the choice of no-slip boundary conditions is too
restrictive to be an accurate description of the physics of chromatin in the
nucleus. After all, the nuclear envelope has a complex nature and is itself
actively fluctuating chu2017origin , which will significantly affect the
possible motions of the fluids within. Furthermore, in this model we are not
allowing any change in the total nuclear volume, or permeation of solvent
through the boundary, both of which occur in real biological systems. These,
and many other improvements on the model, could be added in a systematic way
to this hydrodynamic description of chromatin dynamics, hopefully increasing
our understanding of the complex physics of the genome.

Acknowledgements AZ is grateful for support from the NSF Grants CAREER
PHY-1554880, CMMI-1762506, PHY-2210541, DMS-2153432, and NYU MRSEC
DMR-1420073. This research was supported in part by the National Science
Foundation under Grant No. NSF PHY-1748958. AZ and AYG acknowledge useful
discussions with participants of the 2020 virtual KITP program on ”Biological
Physics of Chromosomes”. Numerical solutions were calculated using New York
University (NYU) High Performance Computing cluster.

Data Availability Data sharing not applicable to this article as no datasets
were generated or analysed during the current study.

##  Appendix A Completeness of solenoidal solutions

In the main text, we simply stated that the solenoidal solutions, that is
those solutions with $\nabla\cdot\mathbf{w}=0$, are spanned by the expansion

$\mathbf{w}_{\perp}=\sum_{lmn}a_{lmn}j_{l}\left(\frac{\alpha_{ln}r}{R}\right)\mathbf{\Phi}_{lm}\;.$
(19)

However, we did not show that there is no combination coming from the basis
functions $\mathbf{\Psi}_{lm},\;\mathbf{Y}_{lm}$. Here we will demonstrate
that this is in fact the case.

Consider one term in the expansion
$\mathbf{w}_{\perp}=E^{r}(r)\mathbf{Y}_{lm}+E^{\theta}(r)\mathbf{\Psi}_{lm}$.
The divergence of this field is

$\begin{split}\nabla\cdot\mathbf{w}_{\perp}&=\left(\partial_{r}E^{r}+\frac{2}{r}E^{r}-\frac{l(l+1)}{r}E^{\theta}\right)Y_{lm}\\\
&=0\;.\end{split}$ (20)

Since the basis functions must be eigenfunctions of the spherical Laplacian
and meet the no-slip boundary condition at $r=R$, we must have

$E^{r}=Cj_{l}\left(\frac{r\alpha_{ln}}{R}\right)$ (21)

for some constant C. Enforcing the divergence-free condition, we solve for
$E^{\theta}$:

$E^{\theta}=C\left(\frac{2+l}{r}j_{l}\left(\frac{\alpha_{ln}r}{R}\right)-\frac{\alpha_{ln}}{R}j_{l+1}\left(\frac{\alpha_{ln}r}{R}\right)\right)$
(22)

which cannot be made to meet the boundary condition at $r=R$. Thus we have
shown that the expansion (19) spans possible divergence-free functions that
meet our boundary conditions and are eigenfunctions of the spherical Laplace
operator.

##  Appendix B Numerical methods

We wrote a numerical scheme in `python` which successfully integrates the
equations of motion

$\begin{split}&\left(1-\lambda_{s}^{2}\nabla\nabla\cdot+\lambda^{2}\nabla\times\nabla\times\right)\tau\partial_{t}\mathbf{w}-2\lambda_{d}^{2}\nabla\nabla\cdot\mathbf{w}\\\
&\quad\quad=\frac{f\rho}{\zeta}\tau\partial_{t}\mathbf{m}\\\
&\tau\partial_{t}\mathbf{m}=2(\mathbf{m}_{eq}(\mathbf{w})-\mathbf{m})\\\
&\mathbf{m}_{eq}=\frac{\mathbf{w}\tau}{3a}\left(1-\frac{3}{5}\left(\frac{\mathbf{w}\tau}{3a}\right)^{2}\right)\;.\end{split}$
(23)

Given the initial fields $w$ and $m$ at time $t=0$, we first compute their
decomposition into normal modes, i.e., eigenfunctions (15). This process is
outlined in further detail below. These modes are uncoupled in the equation of
motion (4) for $\mathbf{w}$, and so we can evolve the normal modes for
$\mathbf{w}$ one time-step directly using equations (16a,16b), discretized
using an Adams-Bashforth two-step scheme. The resulting difference equations
for $a_{nlm},b_{nlm}$ are shown below.

However, at each time-step, the field $\mathbf{m}$ must also be updated,
following equation (8). To do this, the normal modes are reassembled into the
full spatial dependence of $\mathbf{m},\mathbf{w}$, and equation (8) is
evolved directly. Then the procedure above is repeated.

Throughout this scheme, we used the `scipy.special` package for the definition
of the spherical Bessel functions, but the vector spherical harmonics were
implemented through the freely available package `shtns`, detailed in the work
schaeffer2013efficient .

###  B.1 Time evolution

The time evolution was implemented using an Adams-Bashforth two-step scheme
for the source terms, combined with a mixed Crank-Nicolson scheme in the case
of the equation for $b_{nlm}$ to ensure stability (for the details on the
definitions of these schemes, see press1988numerical ):

$\begin{split}&a_{nlm}^{t+1}=\Delta
t\frac{\epsilon+1}{1+\left(\frac{\lambda\alpha_{ln}}{R}\right)^{2}}\left(\frac{3}{2}\left(p_{\perp
nlm}\right)^{t}-\frac{1}{2}\left(p_{\perp nlm}\right)^{t-1}\right)\\\
&b_{nlm}^{t+1}=\Delta t\frac{B_{ln}}{1+\Delta
tA_{ln}}\left(\frac{3}{2}\left(p_{\parallel
nlm}\right)^{t}-\frac{1}{2}\left(p_{\parallel nlm}\right)^{t-1}\right)\\\
&\quad\quad+\frac{1-\Delta tA_{ln}}{1+\Delta tA_{ln}}b_{nlm}^{t-1}\\\
&A_{ln}=\frac{\left(\lambda_{d}\alpha_{ln}/R\right)^{2}}{1+\left(\lambda_{s}\alpha_{ln}/R\right)^{2}},\;B_{ln}=\frac{\epsilon+1}{1+\left(\lambda_{s}\alpha_{ln}/R\right)^{2}}\;,\end{split}$
(24)

where we have defined the field $\mathbf{p}=\partial_{t}\mathbf{m}$ as
shorthand.

###  B.2 Mode decomposition

The remaining task is the conversion from spatial to spherical harmonic
representations, and vice versa. The angular dependence of the decomposition
is taken care of by the routines within the `shtns` package, however we had to
implement our own decomposition of the radial part into Bessel functions.

We will illustrate this whole process with a simple example. Consider a vector
field $\mathbf{f}(\mathbf{r})$, which we seek to decompose into the spherical
harmonic and Bessel function basis, giving coefficients $f_{\perp
nlm},f_{\parallel nlm}$. We perform this decomposition in two steps. First, we
loop over radial slices $\\{r_{i}\\}$, and for each of these values we call
`shtns` to decompose the vector field into spherical components

$\begin{split}\mathbf{f}(r_{i},\theta,\phi)&=\sum_{lm}Q_{lm}(r_{i})\mathbf{Y}_{lm}(\theta,\phi)\\\
&+S_{lm}(r_{i})\mathbf{\Psi}_{lm}(\theta,\phi)\\\
&+T_{lm}(r_{i})\mathbf{\Phi}_{lm}(r_{i})\;,\end{split}$ (25)

which gives us a discrete radial representation of the functions
$Q_{lm}(r),S_{lm}(r),T_{lm}(r)$. We can directly read off that
$T_{lm}(r)=\sum_{nlm}f_{\perp nlm}j_{l}(\alpha_{ln}r/R)$. However, the
definition of the longitudinal modes is more ambiguous. This is resolved by
considering the gradient of a generic scalar field $\rho(\mathbf{r})$ in terms
of VSH:

$\nabla\left(\rho_{lm}(r)Y_{lm}\right)=\frac{\rho_{lm}(r)}{r}\mathbf{\Psi}_{lm}+\rho^{\prime}_{lm}(r)\mathbf{Y}_{lm}$
(26)

Therefore, for a given spherical harmonic mode we must have

$rS_{lm}(r)=\int_{0}^{r}Q_{lm}(r^{\prime})dr^{\prime}=\sum_{nlm}f_{\parallel
nlm}j_{l}(\alpha_{ln}r/R)\;.$ (27)

For the mode where $l=0,m=0$, we decompose $Q_{00}$ to obtain the coefficients
of the derivative of $j_{0}(r\alpha_{0n}/R)$, since $\Psi_{00}=0$. Otherwise
we decompose the function $rS_{lm}(r)$ to get coefficients of
$j_{l}(r\alpha_{ln})$, as this is computationally simpler. Both in the case of
longitudinal and transverse flows, the task is always reduced to the
decomposition of a function defined over the interval $[0,1]$ into spherical
Bessel functions

$\begin{split}f_{l}(x)&=\sum_{n}f_{ln}j_{l}\left(x\alpha_{ln}\right),\;l>0\\\
f_{0}(x)&=\sum_{n}f_{0n}j_{0}^{\prime}\left(n\pi x\right)=-\sum_{n}n\pi
f_{0n}j_{1}(n\pi x)\;.\end{split}$ (28)

where we have suppressed the index $m$ as it does not affect this part of the
calculation, and inserted the known roots $\alpha_{0n}=n\pi$. To find the
components $f_{ln}$, we exploit the following orthogonality relation

$\begin{split}&\int_{0}^{1}j_{l}(x\alpha_{ln})j_{l}(x\alpha_{lm})x^{2}dx=\\\
&\quad\frac{\delta_{nm}}{2}\times\begin{cases}-j_{l-1}(\alpha_{ln})j_{l+1}(\alpha_{ln}),\;l>0\;.\\\
\frac{1}{(n\pi)^{2}},\;l=0\end{cases}\end{split}$ (29)

Multiplying equations 28 by $x^{2}$ and integrating, we obtain

$\begin{split}f_{ln}&=-2\frac{\int_{0}^{1}x^{2}j_{l}(x\alpha_{ln})f_{l}(x)dx}{j_{l-1}(\alpha_{ln})j_{l+1}(\alpha_{ln})},\;l>0\\\
f_{0n}&=-2n\pi\int_{0}^{1}x^{2}j_{1}(n\pi x)f_{0}(x)dx\;.\end{split}$ (30)

These integrals are then implemented in our code. To perform the integrals
over the discrete set of radii $\\{r_{i}\\}$, we used Simpson’s rule as
implemented in `scipy.integrate`.

The inverse operation, going from coefficients to the full spatial dependence,
is far simpler: for each $l$, we sum over Bessel functions multiplied by their
coefficients $f_{ln}$. This gives us the radial dependence of the coefficients
$Q_{lm}(r),S_{lm}(r),T_{lm}(r)$. We then use the `shtns` package to convert
those back to spatial representations, radial shell by radial shell.

## References

  *   * (1) B. Alberts, A. Johnson, J. Lewis, D. Morgan, M. Raff, K. Roberts, P. Walter, _Molecular Biology of the Cell_ (W.W. Norton, 2017)
  * (2) K.E. Van Holde, _Chromatin_ (Springer Science & Business Media, 2012)
  * (3) D.L. Spector, Macromolecular domains within the cell nucleus.  Annu. Rev. Cell Biol. 9(1), 265–315 (1993)
  * (4) K. Rippe, Dynamic organization of the cell nucleus.  Curr. Opin. Genet. Dev. 17(5), 373–380 (2007)
  * (5) A. Zidovska, The rich inner life of the cell nucleus: dynamic organization, active flows, and emergent rheology.  Biophys. Rev. 12(5), 1093–1106 (2020)
  * (6) M.R. Hübner, D.L. Spector, Chromatin dynamics.  Annu. Rev. Biophys. 39, 471–489 (2010)
  * (7) B. Fierz, M.G. Poirier, Biophysics of chromatin dynamics.  Annu. Rev. Biophys. 48, 321–345 (2019)
  * (8) A. Zidovska, The self-stirred genome: large-scale chromatin dynamics, its biophysical origins and implications.  Curr. Opin. Genet. Dev. 61, 83–90 (2020)
  * (9) M.R. Hübner, M.A. Eckersley-Maslin, D.L. Spector, Chromatin organization and transcriptional regulation.  Curr. Opin. Genet. Dev. 23(2), 89–95 (2013)
  * (10) W.A. Bickmore, B. van Steensel, Genome architecture: domain organization of interphase chromosomes.  Cell 152(6), 1270–1284 (2013)
  * (11) L. Mirny, J. Dekker, Mechanisms of chromosome folding and nuclear organization: their interplay and open questions.  Cold Spring Harb. Perspect. Biol. 14(7), a040147 (2022)
  * (12) E. Soutoglou, T. Misteli, Mobility and immobility of chromatin in transcription and genome stability.  Curr. Opin. Genet. Dev. 17(5), 435–442 (2007)
  * (13) C. Lanctôt, T. Cheutin, M. Cremer, G. Cavalli, T. Cremer, Dynamic genome architecture in the nuclear space: regulation of gene expression in three dimensions.  Nat. Rev. Genet. 8(2), 104–115 (2007)
  * (14) J. Miné-Hattab, R. Rothstein, DNA in motion during double-strand break repair.  Trends Cell Biol. 23(11), 529–536 (2013)
  * (15) T. Lebeaupin, H. Sellou, G. Timinszky, S. Huet, Chromatin dynamics at DNA breaks: what, how and why?  AIMS Biophys. 2(4), 458–475 (2015)
  * (16) J.A. Eaton, A. Zidovska, Structural and dynamical signatures of local DNA damage in live cells.  Biophys. J 118(9), 2168–2180 (2020)
  * (17) W. Marshall, A. Straight, J. Marko, J. Swedlow, A. Dernburg, A. Belmont, A. Murray, D. Agard, J. Sedat, Interphase chromosomes undergo constrained diffusional motion in living cells.  Curr. Biol. 7(12), 930–939 (1997)
  * (18) V. Levi, Q. Ruan, M. Plutz, A.S. Belmont, E. Gratton, Chromatin dynamics in interphase cells revealed by tracking in a two-photon excitation microscope.  Biophys. J. 89(6), 4275–4285 (2005)
  * (19) C.H. Chuang, A.E. Carpenter, B. Fuchsova, T. Johnson, P. de Lanerolle, A.S. Belmont, Long-range directional movement of an interphase chromosome site.  Curr. Biol. 16(8), 825–831 (2006)
  * (20) L.A. Mirny, I. Solovei, Keeping chromatin in the loop (s).  Nat. Rev. Molec. Cell Biol. 22(7), 439–440 (2021)
  * (21) I.F. Davidson, J.M. Peters, Genome folding through loop extrusion by SMC complexes.  Nat. Rev. Molec. Cell Biol. 22(7), 445–464 (2021)
  * (22) A. Zidovska, D.A. Weitz, T.J. Mitchison, Micron-scale coherence in interphase chromatin dynamics.  Proc. Natl. Acad. Sci. USA 110(39), 15555–15560 (2013)
  * (23) R. Bruinsma, A.Y. Grosberg, Y. Rabin, A. Zidovska, Chromatin hydrodynamics.  Biophys. J. 106(9), 1871–1881 (2014)
  * (24) D. Saintillan, M.J. Shelley, A. Zidovska, Extensile motor activity drives coherent motions in a model of interphase chromatin.  Proc. Natl. Acad. Sci. USA 115(45), 11442–11447 (2018)
  * (25) A. Mahajan, W. Yan, A. Zidovska, D. Saintillan, M.J. Shelley, Euchromatin activity enhances segregation and compaction of heterochromatin in the cell nucleus.  Phys. Rev. X 12, 041033 (2022)
  * (26) G. Shi, L. Liu, C. Hyeon, D. Thirumalai, Interphase human chromosome exhibits out of equilibrium glassy dynamics.  Nat. Commun. 9(1), 1–13 (2018)
  * (27) L. Liu, G. Shi, D. Thirumalai, C. Hyeon, Chain organization of human interphase chromosome determines the spatiotemporal dynamics of chromatin loci.  PLoS Comput. Biol. 14(12), e1006617 (2018)
  * (28) M. Di Pierro, D.A. Potoyan, P.G. Wolynes, J.N. Onuchic, Anomalous diffusion, spatial coherence, and viscoelasticity from the energy landscape of human chromosomes.  Proc. Natl. Acad. Sci. USA 115(30), 7753–7758 (2018)
  * (29) I. Eshghi, A. Zidovska, A.Y. Grosberg, Symmetry-based classification of forces driving chromatin dynamics.  Soft Matter 18(42), 8134–8146 (2022)
  * (30) I. Eshghi, A. Zidovska, A.Y. Grosberg, Activity-driven phase transition causes coherent flows of chromatin (2023).  arXiv:2305.01052 [cond-mat.soft]
  * (31) M. Doi, A. Onuki, Dynamic coupling between stress and composition in polymer solutions and blends.  J. Phys. II 2(8), 1631–1656 (1992)
  * (32) M.C. Marchetti, J.F. Joanny, S. Ramaswamy, T.B. Liverpool, J. Prost, M. Rao, R.A. Simha, Hydrodynamics of soft active matter.  Rev. Mod. Phys. 85(3), 1143 (2013)
  * (33) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, O. Shochet, Novel type of phase transition in a system of self-driven particles.  Phys. Rev. Lett. 75(6), 1226 (1995)
  * (34) J. Toner, Y. Tu, Long-range order in a two-dimensional dynamical XY model: how birds fly together.  Phys. Rev. Lett. 75(23), 4326 (1995)
  * (35) R.M. Adar, J.F. Joanny, Permeation instabilities in active polar gels.  Phys. Rev. Lett. 127, 188001 (2021)
  * (36) H. Kimura, Y. Tao, R.G. Roeder, P.R. Cook, Quantitation of RNA Polymerase II and its transcription factors in an HeLa cell: little soluble holoenzyme but significant amounts of polymerases attached to the nuclear substructure.  Mol. Cell Biol. 19(8), 5383–5392 (1999)
  * (37) I. Eshghi, J.A. Eaton, A. Zidovska, Interphase chromatin undergoes a local sol-gel transition upon cell differentiation.  Phys. Rev. Lett. 126(22), 228101 (2021)
  * (38) M. Crisp, Q. Liu, K. Roux, J. Rattner, C. Shanahan, B. Burke, P.D. Stahl, D. Hodzic, Coupling of the nucleus and cytoplasm: role of the LINC complex.  J. Cell Biol. 172(1), 41–53 (2006)
  * (39) G.K. Batchelor, _An Introduction to Fluid Dynamics_ (Cambridge University Press, 1967)
  * (40) F.Y. Chu, S.C. Haley, A. Zidovska, On the origin of shape fluctuations of the cell nucleus.  Proc. Natl. Acad. Sci. USA 114(39), 10338–10343 (2017)
  * (41) P. Morse, H. Feshbach, _Methods of Theoretical Physics_ (Feshbach Publishing, 1981)
  * (42) N. Schaeffer, Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations.  Geochem. Geophys. Geosyst. 14(3), 751–758 (2013)
  * (43) W.H. Press, W.T. Vetterling, S.A. Teukolsky, B.P. Flannery, _Numerical Recipes_ (Citeseer, 1988)