# On the density of a set of primes associated to an elliptic curve

Nuno Arala<EMAIL_ADDRESS>
###### Abstract.

For a given elliptic curve $E$ defined over the rationals, we study the
density of primes $p$ satisfying $\mathrm{gcd}(\\#E(\mathbb{F}_{p}),p-1)=1$
and give a conjectural value for this density with strong heuristic evidence
for most elliptic curves, in an appropriate sense.

Keywords: _Elliptic Curves, Rational Points, Galois Representations_

##  1\. Introduction

As motivation for the considerations that follow, we quote the following
result, which gives a sufficient condition for an elliptic curve
$E/\mathbb{Q}$ not to have integral points over any cyclotomic ring
$\mathbb{Z}[\zeta_{p^{n}}]$, where $p$ is a fixed prime and $\zeta_{p^{n}}$
denotes a primitive $p^{n}$-th root of unity.

######  Theorem 1.1.

Let $E/\mathbb{Q}$ be an elliptic curve with $E(\mathbb{Q})=0$ and with
conductor $\mathcal{N}_{E}$. Define

$R_{E}=\\{p\nmid\mathcal{N}_{E}\text{
prime}:\mathrm{gcd}(p(p-1),\\#E(\mathbb{F}_{p}))=1\\}\text{.}$

Then $(E-0)(\mathbb{Z}[\zeta_{p^{n}}])=\emptyset$ for $p\in R_{E}$ and $n\geq
1$.

###### Proof.

This is a particular case of Corollary 4 in [2]. ∎

Here $(E-0)(\mathbb{Z}[\zeta_{p^{n}}])$ denotes the set of points in
$E(\mathbb{Q}(\zeta_{p^{n}}))$ which do _not_ reduce to the identity point $0$
in $E(\mathbb{Z}[\zeta_{p^{n}}]/\mathfrak{P})$ for any prime ideal
$\mathfrak{P}$ of $\mathbb{Z}[\zeta_{p^{n}}]$. If
$y^{2}+a_{1}xy+a_{3}y=x^{3}+a_{2}x^{2}+a_{4}x+a_{6}$ is a minimal Weierstrass
equation for $E/\mathbb{Q}$, then the conclusion of the previous theorem can
be phrased by stating that the equation
$y^{2}+a_{1}xy+a_{3}y=x^{3}+a_{2}x^{2}+a_{4}x+a_{6}$ has no solutions in
$(x,y)\in(\mathbb{Z}[\zeta_{p^{n}}])^{2}$.

The previous theorem makes it natural to attempt to estimate the density of
primes contained in $R_{E}$, for a fixed $E$. Obviously the condition
$p\nmid\mathcal{N}_{E}$ is unimportant for density considerations since it
only affects a finite set of primes. Moreover the set of primes such $p$ such
that $\mathrm{gcd}(p,\\#E(\mathbb{F}_{p}))>1$ (i.e. such that
$p\mid\\#E(\mathbb{F}_{p})$) has density $0$. This can be seen as follows. It
is well-known that

$\\#E(\mathbb{F}_{p})=p+1-a_{p}$

where $a_{p}$ is the trace of the Frobenius map acting on the reduction of $E$
over $\mathbb{F}_{p}$. The Hasse bound ([1], Theorem V.1.1) states that
$|a_{p}|<2\sqrt{p}$. On the other hand, if $p$ divides $\\#E(\mathbb{F}_{p})$,
then we must have $a_{p}\equiv 1\pmod{p}$, and for $p\geq 7$ this together
with $|a_{p}|<2\sqrt{p}$ forces $a_{p}=1$. Now the Sato-Tate Conjecture ([5],
[6], [7]) easily implies that the set of primes $p$ for which $a_{p}=1$ has
density $0$. We reduce therefore our question to the following problem.

Problem. Given an elliptic curve $E/\mathbb{Q}$, evaluate the density of
primes $p$ satisfying

$\mathrm{gcd}(p-1,\\#E(\mathbb{F}_{p}))=1\text{.}$

In what follows, we give strong evidence (namely, Theorem 4.4, which is the
main result in the paper) towards the following conjectural (partial) answer
to our problem (for the definition of a Serre curve, see Proposition 4.3 and
the discussion after that).

######  Conjecture 1.2.

Define

$C=\prod_{\ell\text{
prime}}\left(1-\frac{\ell}{(\ell-1)^{2}(\ell+1)}\right)\text{.}$

Let $E/\mathbb{Q}$ be a Serre curve, and let $D$ be the discriminant of the
number field $\mathbb{Q}(\sqrt{\Delta})$, where $\Delta$ is the discriminant
of $E$. Then the density of primes $p$ for which
$\mathrm{gcd}(\\#E(\mathbb{F}_{p}),p-1)=1$ is

$\begin{cases}C&\text{ if }D\equiv 0\pmod{4}\\\ \left(1+\prod_{\ell\mid
D}\frac{-\ell}{\ell^{3}-\ell^{2}-2\ell+1}\right)C&\text{ if }D\equiv
1\pmod{4}\text{.}\end{cases}$

##  2\. Notation

Our approach will make use of Galois representations associated to the torsion
of $E$. We fix the corresponding notation here. For a (not necessarily
positive) integer $N$, we use the standard notation $E[N]$ for the subgroup of
$N$-torsion points of $E$; then $E[N]\cong(\mathbb{Z}/N\mathbb{Z})^{2}$ ([1],
Corollary III.6.4). The absolute Galois group
$G_{\mathbb{Q}}={\mathrm{Gal}}(\overline{\mathbb{Q}}/\mathbb{Q})$ acts on
$E[N]$ via $\mathbb{Z}/N\mathbb{Z}$-automorphisms, yielding a group
homomorphism

$\rho_{E,N}:G_{\mathbb{Q}}\to\mathrm{GL}_{2}(\mathbb{Z}/N\mathbb{Z})\text{.}$

Similarly, denoting by $\mathbb{Q}(E[N])$ the smallest extension of
$\mathbb{Q}$ over which all points of $E[N]$ are rational, the Galois group
${\mathrm{Gal}}(\mathbb{Q}(E[N])/\mathbb{Q})$ acts on $E[N]$ via
$\mathbb{Z}/N\mathbb{Z}$-automorphisms; this yields a group homomorphism
$\overline{\rho}_{E,N}:{\mathrm{Gal}}(\mathbb{Q}(E[N])/\mathbb{Q})\to\mathrm{GL}_{2}(\mathbb{Z}/N\mathbb{Z})$.
These two homomorphisms fit in the following commutative diagram:

${G_{\mathbb{Q}}}$${{\mathrm{Gal}(\mathbb{Q}(E[N])/\mathbb{Q})}}$${\mathrm{GL}_{2}(\mathbb{Z}/N\mathbb{Z})}$$\scriptstyle{\rho_{E,N}}$$\scriptstyle{\overline{\rho}_{E,N}}$

By taking the inverse limit, the maps $(\rho_{E,N})_{N>0}$ yield a
homomorphism

$\rho_{E}:G_{\mathbb{Q}}\to\varprojlim\mathrm{GL}_{2}(\mathbb{Z}/N\mathbb{Z})=\mathrm{GL}_{2}(\widehat{\mathbb{Z}})\text{.}$

##  3\. The local count

The condition $\mathrm{gcd}(p-1,\\#E(\mathbb{F}_{p}))=1$ can be restated as
follows: there is no prime $\ell$ satisfying $p\equiv 1\pmod{\ell}$ and
$\ell\mid\\#E(\mathbb{F}_{p})$. Therefore we begin by evaluating, for a fixed
prime $\ell$, the density of primes $p$ satisfying $p\equiv 1\pmod{\ell}$ and
$\ell\mid\\#E(\mathbb{F}_{p})$.

Recall the Galois representation $\overline{\rho}_{E,\ell}$, which embeds
${\mathrm{Gal}}(\mathbb{Q}(E[\ell])/\mathbb{Q})$ into
$\mathrm{GL}_{2}(\mathbb{F}_{\ell})$; we will work now under the assumption
that this embedding is actually an isomorphism. The main result in this
section is the following.

######  Lemma 3.1.

Let $E/\mathbb{Q}$ be an elliptic curve and let $\ell$ be a prime. Suppose
that $\overline{\rho}_{E,\ell}$ (or, equivalently, $\rho_{E,\ell}$) is
surjective. Then the density of primes $p$ for which $\ell$ is not a common
factor of $p-1$ and $\\#E(\mathbb{F}_{p})$ is

$1-\frac{\ell}{(\ell-1)^{2}(\ell+1)}\text{.}$

In order to prove this we need two auxiliary propositions.

######  Proposition 3.2.

Let $E/\mathbb{Q}$ be an elliptic curve and let $\ell$ be a prime. Let
$p\neq\ell$ be a prime of good reduction for $E$, let $\mathfrak{P}$ denote a
prime above $p$ in $\mathbb{Q}(E[\ell])$ and let $\sigma_{\mathfrak{P}}$
denote the corresponding Frobenius element. Then $\ell$ divides
$\\#E(\mathbb{F}_{p})$ if and only if $1$ is an eigenvalue of
$\overline{\rho}_{E,\ell}(\sigma_{\mathfrak{P}})$.

###### Proof.

Assume first that $1$ is an eigenvalue of
$\overline{\rho}_{E,\ell}(\sigma_{\mathfrak{P}})$. Then
$\sigma_{\mathfrak{P}}$ has a fixed point, say $Q$, in $E[\ell]$. Since
$\sigma_{\mathfrak{P}}(Q)=Q$, reducing this modulo $\mathfrak{P}$ yields
$\phi_{p}(Q^{\prime})=Q^{\prime}$, where $Q^{\prime}$ is the reduction of $Q$
modulo $\mathfrak{P}$ and $\phi_{p}$ is the Frobenius automorphism of
$\overline{\mathbb{F}_{p}}$ given by raising to the $p$-th power. This
equality implies that $Q^{\prime}\in E(\mathbb{F}_{p})$. Since $Q\in E[\ell]$
and the reduction map restricted to $E[\ell]$ is injective ([1], Proposition
IV.3.1) it follows that $Q^{\prime}$ has order $\ell$. Therefore
$E(\mathbb{F}_{p})$ has a point of order $\ell$. Hence by Lagrange’s Theorem
$\ell\mid\\#E(\mathbb{F}_{p})$.

Conversely, assume that $\ell\mid\\#E(\mathbb{F}_{p})$. Then by Cauchy’s
Theorem $E(\mathbb{F}_{p})$ contains a point $Q^{\prime}$ of order $\ell$. We
lift $Q^{\prime}$ to a point $Q$ in $E[\ell]$. Then the fact that
$\phi_{p}(Q^{\prime})=Q^{\prime}$ together with the definition of the
Frobenius element $\sigma_{\mathfrak{P}}$ implies

$\sigma_{\mathfrak{P}}(Q)\equiv Q\pmod{\mathfrak{P}}\text{.}$

But injectivity of reduction modulo $\mathfrak{P}$ on $E[\ell]$ implies that
we actually have $\sigma_{\mathfrak{P}}(Q)=Q$. Therefore
$\sigma_{\mathfrak{P}}$ has a fixed point in $E[\ell]$, which means precisely
that $\overline{\rho}_{E,\ell}(\sigma_{\mathfrak{P}})$ has $1$ as an
eigenvalue. ∎

Similarly, we now express the condition that $\ell\mid p-1$ in terms of the
image of $\sigma_{\mathfrak{P}}$ under the Galois representation
$\overline{\rho}_{E,\ell}$.

######  Proposition 3.3.

Let $E/\mathbb{Q}$ be an elliptic curve and let $\ell$ be a prime. Let
$p\neq\ell$ be a prime, let $\mathfrak{P}$ denote a prime above $p$ in
$\mathbb{Q}(E[\ell])$ and let $\sigma_{\mathfrak{P}}$ denote the corresponding
Frobenius element. Then $\ell$ divides $p-1$ if and only if
$\mathrm{det}(\overline{\rho}_{E,\ell}(\sigma_{\mathfrak{P}}))=1$ (i.e. if
$\overline{\rho}_{E,\ell}(\sigma_{\mathfrak{P}})\in\mathrm{SL}_{2}(\mathbb{F}_{\ell})$).

###### Proof.

We first recall that $\mathbb{Q}(\zeta_{\ell})\subseteq\mathbb{Q}(E[\ell])$,
where $\zeta_{\ell}$ denotes a primitive $\ell$-th root of unity. This is a
consequence of basic properites of the Weil pairing (see [1], section III.8).
Therefore $\mathrm{Gal}(\mathbb{Q}(E[\ell])/\mathbb{Q})$ surjects naturally
onto
$\mathrm{Gal}(\mathbb{Q}(\zeta_{\ell})/\mathbb{Q})\cong\mathbb{F}_{\ell}^{\times}$.
Under $\overline{\rho}_{E,\ell}$,
${\mathrm{Gal}}(\mathbb{Q}(E[\ell])/\mathbb{Q})$ is identified with a subgroup
of $\mathrm{GL}_{2}(\mathbb{F}_{\ell})$, so we get a map from a subgroup of
$\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ to $\mathbb{F}_{\ell}^{\times}$. It is a
standard fact that this map is given by the determinant. In other words, we
have a commutative diagram:

${{\mathrm{Gal}(\mathbb{Q}(E[\ell])/\mathbb{Q})}}$${\mathrm{GL}_{2}(\mathbb{F}_{\ell})}$${\mathbb{F}_{\ell}^{\times}}$$\scriptstyle{\overline{\rho}_{E,\ell}}$$\scriptstyle{\mathrm{det}}$

This reduces our task to proving that $p\equiv 1\pmod{\ell}$ if and only if
$\sigma_{\mathfrak{P}}(\zeta_{\ell})=\zeta_{\ell}$. Since
$\sigma_{\mathfrak{P}}(\zeta_{\ell})=\zeta_{\ell}^{p}$, this is clear. ∎

We are now in a position to prove Lemma 3.1.

###### Proof of Lemma 3.1.

Keeping the notation used in Lemmas 3.2 and 3.3, combining both these lemmas
we see that we are interested in evaluating the density of primes $p$ for
which $\overline{\rho}_{E,\ell}(\sigma_{\mathfrak{P}})$ has determinant $1$
and has $1$ as an eigenvalue (this is the density of the set of “bad” primes,
those for which $\ell$ is a common factor of $p-1$ and
$\\#E(\mathbb{F}_{p})$). In other words, we want primes $p$ such that
$\overline{\rho}_{E,\ell}(\sigma_{\mathfrak{P}})$ is conjugate to a matrix of
the form

$\begin{pmatrix}1&\ast\\\ 0&1\end{pmatrix}\text{.}$

We count these matrices: apart from the identity matrix, all these matrices
fall into the conjugacy class of $T=\begin{pmatrix}1&1\\\ 0&1\end{pmatrix}$.
In fact, for $a\in\mathbb{F}_{\ell}^{\star}$, a calculation shows that

$\begin{pmatrix}1&0\\\ 0&a\end{pmatrix}^{-1}\begin{pmatrix}1&1\\\
0&1\end{pmatrix}\begin{pmatrix}1&0\\\ 0&a\end{pmatrix}=\begin{pmatrix}1&a\\\
0&1\end{pmatrix}\text{.}$

We are then interested in computing the cardinalities of the conjucacy classes
in $\mathbb{GL}_{2}(\mathbb{F}_{\ell})$ of the identity and $T$. The former
has cardinality $1$; we now compute the order of the stabilizer of $T$ under
the conjugation action. A matrix $M=\begin{pmatrix}a&c\\\
b&d\end{pmatrix}\in\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ stabilizes $T$ if and
only if

$\begin{pmatrix}a&c\\\ b&d\end{pmatrix}\begin{pmatrix}1&1\\\
0&1\end{pmatrix}=\begin{pmatrix}1&1\\\ 0&1\end{pmatrix}\begin{pmatrix}a&c\\\
b&d\end{pmatrix}\text{,}\quad\text{i.e.,}\quad\begin{pmatrix}a&a+c\\\
b&b+d\end{pmatrix}=\begin{pmatrix}a+b&c+d\\\ b&d\end{pmatrix}\text{.}$

The previous equality is equivalent to $a=d$ and $b=0$. This gives $\ell^{2}$
possible matrices, but those $\ell$ of them for which $a=d=0$ are not in
$\mathrm{GL}_{2}(\mathbb{F}_{\ell})$, so we conclude that the stabilizer of
$T$ has $\ell^{2}-\ell$ elements. By the orbit-stabilizer theorem together
with the well-known fact that
$\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})=(\ell^{2}-1)(\ell^{2}-\ell)$, we
conclude that the conjugacy class of $T$ has

$\frac{\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})}{\ell^{2}-\ell}=\ell^{2}-1$

elements. Therefore there are $\ell^{2}$ matrices in
$\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ which have determinant $1$ and have $1$
as eigenvalue.

Now we use Chebotarev’s density theorem: by assumption
$\overline{\rho}_{E,\ell}$ identifies
${\mathrm{Gal}}(\mathbb{Q}(E[\ell])/\mathbb{Q})$ with
$\mathrm{GL}_{2}(\mathbb{F}_{\ell})$, and we want to count primes $p$ for
which $\overline{\rho}_{E,\ell}(\sigma_{\mathfrak{P}})$ falls outside one of
two conjugacy classes of $\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ under this
identification. By the computation above, these conjugacy classes altogether
have $\ell^{2}$ elements, and since $\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ has
$(\ell^{2}-1)(\ell^{2}-\ell)$ elements, the remaining conjugacy classes
comprise $(\ell^{2}-1)(\ell^{2}-\ell)-\ell^{2}$ elements of
$\mathrm{GL}_{2}(\mathbb{F}_{\ell})$. By Chebotarev’s density theorem it
follows that the desired density is

$\frac{(\ell^{2}-1)(\ell^{2}-\ell)-\ell^{2}}{(\ell^{2}-1)(\ell^{2}-\ell)}=1-\frac{\ell}{(\ell-1)^{2}(\ell+1)}\text{,}$

as we wanted. ∎

##  4\. The global count

Recall that, given an elliptic curve $E/\mathbb{Q}$, we are interested in
computing the density of rational primes $p$ satisfying
$\mathrm{gcd}(p-1,\\#E(\mathbb{F}_{p}))=1$. A prime $p$ satisfies this if and
only if, for every prime $\ell$, the conditions $p\equiv 1\pmod{\ell}$ and
$\ell\mid\\#E(\mathbb{F}_{p})$ do not hold simultaneously. Suppose that the
Galois representations $\overline{\rho}_{E,\ell}$ are surjective for all
primes $\ell$. A naive probabilistic reasoning, based on Lemma 3.1, then
suggests that the desired density is

$\prod_{\ell\text{
prime}}\left(1-\frac{\ell}{(\ell-1)^{2}(\ell+1)}\right)\approx
0.24238005\text{.}$ (4.1)

This naive reasoning, however, turns out to give the wrong answer, and the
reason is that it assumes some sort of independence between the events

$(p\equiv 1\pmod{\ell}\text{ and }\ell\mid\\#E(\mathbb{F}_{p}))$

as $\ell$ ranges over the primes. This does not necessarily hold, and we need
to add some factors to our prediction (4.1) (the so-called entanglement
correction factors) to account for this lack of independence.

To explain this lack of independence, let $\ell_{1}$ and $\ell_{2}$ be
distinct primes. As we saw before, in Lemmas 3.2 and 3.3, whether $\ell_{1}$
divides $\\#E(\mathbb{F}_{p})$ is determined by the Frobenius over $p$ in
${\mathrm{Gal}}(\mathbb{Q}(E[\ell_{1}])/\mathbb{Q})$, and, similarly, whether
$\ell_{2}$ divides $\\#E(\mathbb{F}_{p})$ is determined by the Frobenius over
$p$ in ${\mathrm{Gal}}(\mathbb{Q}(E[\ell_{2}])/\mathbb{Q})$. Denote these two
automorphisms by $\sigma_{p,\ell_{1}}$ and $\sigma_{p,\ell_{2}}$,
respectively. Then $\sigma_{p,\ell_{1}}$ and $\sigma_{p,\ell_{2}}$ are the
restrictions to $\mathbb{Q}(E[\ell_{1}])$ and $\mathbb{Q}(E[\ell_{2}])$,
respectively, of an element of
${\mathrm{Gal}}(\mathbb{Q}(E[\ell_{1}\ell_{2}])/\mathbb{Q})$, namely the
Frobenius $\sigma_{p,\ell_{1}\ell_{2}}$ over $p$ in
${\mathrm{Gal}}(\mathbb{Q}(E[\ell_{1}\ell_{2}])/\mathbb{Q})$.

${{\mathrm{Gal}(\mathbb{Q}(E[\ell_{1}\ell_{2}])/\mathbb{Q})}}$${{\sigma_{p,\ell_{1}\ell_{2}}}}$${{\mathrm{Gal}(\mathbb{Q}(E[\ell_{1}])/\mathbb{Q})\quad}}$${{\quad\mathrm{Gal}(\mathbb{Q}(E[\ell_{2}])/\mathbb{Q})}}$${{\sigma_{p,\ell_{1}}}}$${{\sigma_{p,\ell_{2}}}}$${{\mathrm{Gal}(\mathbb{Q}(E[\ell_{1}])\cap\mathbb{Q}(E[\ell_{2}])/\mathbb{Q})}}$${\sigma}$

But this implies that $\sigma_{p,\ell_{1}}$ and $\sigma_{p,\ell_{2}}$ restrict
to the same element of
$\mathrm{Gal}(\mathbb{Q}(E[\ell_{1}])\cap\mathbb{Q}(E[\ell_{2}])/\mathbb{Q})$.
If the intersection field $\mathbb{Q}(E[\ell_{1}])\cap\mathbb{Q}(E[\ell_{2}])$
is strictly larger than $\mathbb{Q}$, then this gives a non-trivial dependence
between $\sigma_{p,\ell_{1}}$ and $\sigma_{p,\ell_{2}}$, for an arbitrary
prime $p$.

######  Example 4.1.

Suppose $E/\mathbb{Q}$ is an elliptic curve whose discriminant $\Delta$ is of
the form $-3n^{2}$, for some integer $n$, and that $\overline{\rho}_{E,2}$ and
$\overline{\rho}_{E,3}$ are surjective. Then, according to the naive heuristic
that predicted (4.1), the density of primes $p$ for which $2$ and $3$ do not
divide $\\#E(\mathbb{F}_{p})$ should be

$\prod_{\ell=2,3}\left(1-\frac{\ell}{(\ell-1)^{2}(\ell+1)}\right)=\frac{13}{48}\text{.}$
(4.2)

However, this is not correct, because the fields $\mathbb{Q}(E[2])$ and
$\mathbb{Q}(E[3])$ intersect non-trivially: on the one hand it is well-known
that $\mathbb{Q}(E[2])$ contains $\sqrt{\Delta}$ (and hence contains
$\sqrt{-3}$, given our assumption on $\Delta$), and on the other hand, as
remarked before, $\mathbb{Q}(E[3])$ contains $\zeta_{3}=\frac{1+\sqrt{-3}}{2}$
(and hence contains $\sqrt{-3}$). Therefore $\mathbb{Q}(E[2])$ and
$\mathbb{Q}(E[3])$ both contain $\mathbb{Q}(\sqrt{-3})$:

${{\mathbb{Q}(E[6])}}$${{\mathbb{Q}(E[2])}}$${{\mathbb{Q}(E[3])}}$${\mathbb{Q}(\sqrt{-3})}$${\mathbb{Q}}$

This puts a restriction on the possible pairs $(\sigma_{p,2},\sigma_{p,3})$,
since for any prime $p$ both of these must act on $\sqrt{-3}$ in the same way.
The computations that follow imply that (if $\mathbb{Q}(\sqrt{-3})$ is the
full intersection $\mathbb{Q}(E[2])\cap\mathbb{Q}(E[3])$) the correct constant
in place of (4.2) is $\frac{5}{24}$.

More generally, if $\ell_{1},\ldots,\ell_{m}$ and
$\tilde{\ell}_{1},\ldots,\tilde{\ell}_{n}$ are distinct primes, then if

$\mathbb{Q}(E[\ell_{1}\cdots\ell_{m}])\cap\mathbb{Q}(E[\tilde{\ell}_{1}\cdots\tilde{\ell}_{n}])$
(4.3)

is strictly larger than $\mathbb{Q}$, then the $m$-tuple of Frobenius
automorphisms $(\sigma_{p,\ell_{1}},\ldots,\sigma_{p,\ell_{m}})$ places a non-
trivial restriction on the $n$-tuple
$(\sigma_{p,\tilde{\ell}_{1}},\ldots,\sigma_{p,\tilde{\ell}_{n}})$. In
general, the correction factors we need to multiply (4.1) with will thus
depend on how intersections of the form (4.3) behave as
$\ell_{1},\ldots,\ell_{m}$ and $\tilde{\ell}_{1},\ldots,\tilde{\ell}_{n}$
range over the primes. On what follows we will give a corrected version of
(4.1) for so-called Serre curves, the definition of which we give below.
Recall the homomorphism

$\rho_{E}:G_{\mathbb{Q}}\to\varprojlim\mathrm{GL}_{2}(\mathbb{Z}/N\mathbb{Z})=\mathrm{GL}_{2}(\widehat{\mathbb{Z}})$

(see Section 2). This homomorphism is _never_ surjective, as a theorem of
Serre (see [3]) shows; indeed we show below, following Serre, that there is an
index $2$ subgroup of $\mathrm{GL}_{2}(\widehat{\mathbb{Z}})$ in which
$\rho_{E}(G_{\mathbb{Q}})$ is contained. The failure of $\rho_{E}$ to surject
is in some cases explained by non-trivial intersections of the form (4.3), as
we see in the proof below. We introduce a convenient piece of notation first.

######  Definition 4.2.

Given $\gamma\in\mathrm{GL}_{2}(\mathbb{F}_{2})$, let it act on
$\mathbb{F}_{2}^{2}$ in the natural way; it then permutes the three nonzero
elements of $\mathbb{F}_{2}^{2}$. We define ${\mathrm{sgn}}(\gamma)$ to be the
sign of this permutation.

######  Proposition 4.3.

Let $E/\mathbb{Q}$ be an elliptic curve, let $\Delta$ be its discriminant, and
let $D$ be the discriminant of the $\mathbb{Q}(\sqrt{\Delta})$. Then,

  1. (i)

If $D\equiv 0\pmod{4}$, then there exists a primitive character $\chi$ modulo
$D$ such that

$\rho_{E,D}(G_{\mathbb{Q}})\subseteq\\{M\in\mathrm{GL}_{2}(\mathbb{Z}/D\mathbb{Z}):\chi(\mathrm{det}(M))={\mathrm{sgn}}(M\bmod{2})\\}\text{,}$

where the right hand side has index $2$ in
$\mathrm{GL}_{2}(\mathbb{Z}/D\mathbb{Z})$.

  2. (ii)

If $D\equiv 1\pmod{4}$, then there exists a primitive character $\chi$ modulo
$D$ such that

$\rho_{E,2D}(G_{\mathbb{Q}})\subseteq\\{M\in\mathrm{GL}_{2}(\mathbb{Z}/2D\mathbb{Z}):\chi(\mathrm{det}(M))={\mathrm{sgn}}(M\bmod{2})\\}\text{,}$

where the right hand side has index $2$ in
$\mathrm{GL}_{2}(\mathbb{Z}/2D\mathbb{Z})$.

Therefore, in either case, there is a primitive character $\chi$ modulo $D$
such that

$\rho_{E}(G_{\mathbb{Q}})\subseteq G_{E}$ (4.4)

where

$G_{E}=\\{M\in\mathrm{GL}_{2}(\widehat{\mathbb{Z}}):\chi(\mathrm{det}(M\pmod{D}))={\mathrm{sgn}}(M\bmod{2})\\}$

has index $2$ in $\mathrm{GL}_{2}(\widehat{\mathbb{Z}})$.

###### Proof.

We prove only the case $D\equiv 1\pmod{4}$, because it is the one that matters
the most to us and the proof in the case $D\equiv 0\pmod{4}$ is completely
analogous. We consider the following diagram of field extensions:

${{\mathbb{Q}(E[2D])}}$${{\mathbb{Q}(E[D])}}$${{\mathbb{Q}(E[2])}}$${\mathbb{Q}(\zeta_{D})}$${\mathbb{Q}(\sqrt{D})}$${\mathbb{Q}}$

Here we used the standard inclusion
$\mathbb{Q}(\sqrt{D})\subseteq\mathbb{Q}(\zeta_{D})$, that holds for any
quadratic field discriminant $D$. Take now any $\sigma\in G_{\mathbb{Q}}$, and
let $M=\rho_{E,2D}(\sigma)$. We now describe the action of $\sigma$ on
$\sqrt{D}$ in terms of $M$:

  * •

On the one hand, if $y^{2}=x^{3}+Ax+B$ is a reduced Weierstrass equation for
$E$, then the $2$-torsion points of $E$ are $0$, $(e_{1},0)$, $(e_{2},0)$ and
$(e_{3},0)$ where $e_{1},e_{2},e_{3}$ are the roots of $x^{3}+Ax+B$. Since
$\sqrt{D}$ is $(e_{1}-e_{2})(e_{1}-e_{3})(e_{2}-e_{3})$ up to a rational
factor, $\sigma(\sqrt{D})$ equals $\sqrt{D}$ or $-\sqrt{D}$ according to
whether $\sigma$ permutes $e_{1},e_{2},e_{3}$ via an even permutation or an
odd permutation. The sign of this permutation is what we defined as the sign
of $M\bmod{2}$:

$\sigma(\sqrt{D})=\begin{cases}\sqrt{D}&\text{ if
}{\mathrm{sgn}}(M\bmod{2})=1\\\ -\sqrt{D}&\text{ if
}{\mathrm{sgn}}(M\bmod{2})=-1\text{.}\end{cases}$

  * •

On the other hand, as we showed in the proof of 3.3, if we identify
${\mathrm{Gal}}(\mathbb{Q}(E[D])/\mathbb{Q})$ with
$\mathrm{GL}_{2}(\mathbb{Z}/D\mathbb{Z})$ in the natural way then $\sigma$
acts on $\mathbb{Q}(\zeta_{D})$ via the automorphism corresponding to
$\mathrm{det}(M)$ (we showed this specifically for prime $D$, but the proof
carries through in general). We now identify
${\mathrm{Gal}}(\mathbb{Q}(\sqrt{D})/\mathbb{Q})$ with $\\{\pm 1\\}$ and
consider the map $\chi$ fitting in the following commutative square:

${\mathrm{Gal}(\mathbb{Q}(\zeta_{D})/\mathbb{Q})}$${\mathrm{Gal}(\mathbb{Q}(\sqrt{D})/\mathbb{Q})}$${(\mathbb{Z}/D\mathbb{Z})^{\times}}$${\\{\pm
1\\}}$$\scriptstyle{\cong}$$\scriptstyle{\cong}$$\scriptstyle{\chi}$

We can then regard $\chi$ as a Dirichlet character modulo $D$, and its
primitivity is an easy consequence of the fact that $|D|$ is actually the
_smallest_ positive integer $m$ such that
$\mathbb{Q}(\sqrt{D})\subseteq\mathbb{Q}(\zeta_{m})$. Now $\sigma$ fixes or
moves $\sqrt{D}$ according to whether $\chi(\mathrm{det}(M))$ is $1$ or $-1$
respectively, so that

$\sigma(\sqrt{D})=\begin{cases}\sqrt{D}&\text{ if }\chi(\mathrm{det}(M))=1\\\
-\sqrt{D}&\text{ if }\chi(\mathrm{det}(M))=-1\text{.}\end{cases}$

By comparing the two points above it follows immediately that

${\mathrm{sgn}}(M\bmod{2})=\chi(\mathrm{det}(M))\text{,}$

as desired. It is clear that exactly half of the elements of
$\mathrm{GL}_{2}(\mathbb{Z}/D\mathbb{Z})$ have this property. ∎

We say that an elliptic curve $E/\mathbb{Q}$ is a _Serre curve_ if equality
holds in (4.4). As it turns out, this condition holds quite often; in [4], N.
Jones proves that almost all elliptic curves over $\mathbb{Q}$ are Serre
curves (for a precise statement see Theorem 4 in [4]). This makes it
reasonable to focus our attention henceforth on Serre curves, which are,
according to Proposition 4.3, those for which $\rho_{E}(G_{\mathbb{Q}})$ is as
large as possible. Note that for a Serre curve the maps
$\overline{\rho}_{E,\ell}$, for prime $\ell$, are automatically surjective.

######  Theorem 4.4.

Let $E/\mathbb{Q}$ be a Serre curve, and keep the notation from Proposition
4.3. Let $S$ be a finite set of primes, and define

$\delta_{S}=\prod_{\ell\in
S}\left(1-\frac{\ell}{(\ell-1)^{2}(\ell+1)}\right)\text{.}$

Then,

  1. (i)

If $D\equiv 0\pmod{4}$, the density of primes $p$ for which the conditions
$\ell\mid\\#E(\mathbb{F}_{p})$ and $p\equiv 1\pmod{\ell}$ do not hold
simultaneously for any $\ell\in S$ is $\delta_{S}$.

  2. (ii)

If $D\equiv 1\pmod{4}$, the density of primes $p$ for which the conditions
$\ell\mid\\#E(\mathbb{F}_{p})$ and $p\equiv 1\pmod{\ell}$ do not hold
simultaenously for any $\ell\in S$ is

$\begin{cases}\delta_{S}&\text{ if }2D\nmid\prod_{\ell\in S}\ell\\\
\left(1+\prod_{\ell\mid
D}\frac{-\ell}{\ell^{3}-\ell^{2}-2\ell+1}\right)\delta_{S}&\text{ if
}2D\mid\prod_{\ell\in S}\ell\text{.}\end{cases}$

###### Proof.

Denote by $\Pi$ the product $\prod_{\ell\in S}\ell$. If the natural composite
injection

${{\mathrm{Gal}(\mathbb{Q}(E[\Pi])/\mathbb{Q})}}$${{\prod_{\ell\in
S}\mathrm{Gal}(\mathbb{Q}(E[\ell])/\mathbb{Q})}}$${\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell})}$$\scriptstyle{\cong}$

is an isomorphism, then Lemma 3.1 extends multiplicatively: by Propositions
3.2 and 3.3 and Chebotarev’s density theorem, we are reduced to counting
tuples consisting of one element of $\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ for
each $\ell\in S$, such that the given element of
$\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ is not conjugate to a matrix of the form

$\begin{pmatrix}1&\ast\\\ 0&1\end{pmatrix}\text{.}$

Thus the required density is the product of local densities as computed in the
proof of Lemma 3.1, i.e. $\delta_{S}$. We will then show that the map above is
an isomorphism when $D\equiv 0\pmod{4}$ and also when $D\equiv 1\pmod{4}$ and
not all prime factors of $2D$ are in $S$, establishing part (i) and the first
case of part (ii).

Suppose $D\equiv 0\pmod{4}$; let $2^{k}$ be the highest power of $2$ that
divides $D$. Then $k\geq 2$ (in fact $k\in\\{2,3\\}$). Then the Dirichlet
character $\chi$ from Proposition 4.3 factors as

$\chi_{2^{k}}\prod_{\begin{subarray}{c}\ell\mid D\\\ \ell\neq
2\end{subarray}}\chi_{\ell}$

where $\chi_{2^{k}}$ is a real primitive Dirichlet character modulo $2^{k}$
and $\chi_{\ell}$ is a real primitive Dirichlet character modulo $\ell$ for
each odd prime $\ell$ dividing $D$. Take a tuple

$(M_{\ell})_{\ell\in S}\in\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell})\text{;}$

we want to show that this tuple lies in the image of $\rho_{E,\Pi}$, where we
identified $\mathrm{GL}_{2}(\Pi)$ with $\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell})$. Extend this tuple to a tuple
$M=(M_{\ell})_{\ell\in S\cup T}\in\prod_{\ell\in S\cup
T}\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ in an arbitrary way; it suffices to show
that this tuple lies in the image of $\rho_{E,\Pi^{\prime}}$, where
$\Pi^{\prime}$ is the product of the primes in $S\cup T$. According to the
Serre condition in Proposition 4.3, it suffices to show that there exists a
lift $M_{2^{k}}$ of $M_{2}$ in $\mathrm{GL}_{2}(\mathbb{Z}/2^{k}\mathbb{Z})$
such that

$\chi_{2^{k}}(\mathrm{det}(M_{2^{k}}))\prod_{\ell\in
T\setminus\\{2\\}}\chi_{\ell}(\mathrm{det}(M_{\ell}))={\mathrm{sgn}}(M_{2})\text{.}$
(4.5)

Take any lift

$M_{2^{k}}=\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\mathrm{GL}_{2}(\mathbb{Z}/2^{k}\mathbb{Z})\text{.}$

If (4.5) holds we are done. Otherwise take
$\omega\in(\mathbb{Z}/2^{k}\mathbb{Z})^{\times}$ such that
$\chi_{2^{k}}(\omega)=-1$, and take

$M^{\prime}_{2^{k}}=\begin{pmatrix}a&\omega b\\\ c&\omega
d\end{pmatrix}\text{.}$

Then $M^{\prime}_{2^{k}}$ is another lift of $M_{2}$ but
$\mathrm{det}(M^{\prime}_{2^{k}})=\omega\mathrm{det}(M_{2^{k}})$, so
$\chi_{2^{k}}(\mathrm{det}(M^{\prime}_{2^{k}}))=-\chi_{2^{k}}(\mathrm{det}(M_{2^{k}}))$,
and (4.5) holds with $M_{2^{k}}$ replaced by $M^{\prime}_{2^{k}}$.

Suppose now $D\equiv 1\pmod{4}$ and not all prime factors of $2D$ are in $S$.
This time the character $\chi$ factors as

$\chi=\prod_{\ell\mid D}\chi_{\ell}$

where each $\chi_{\ell}$ is a real primitive character modulo $\ell$ for each
$\ell\mid D$. Once more, let

$(M_{\ell})_{\ell\in S}\in\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell})\text{;}$

be a tuple we want to show is in the image of $\rho_{E,\Pi}$, where again we
identified $\mathrm{GL}_{2}(\Pi)$ with $\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ via the Chinese Remainder Theorem. Let
now $T$ be the set of prime factors of $2D$, and let $q$ be a prime factor of
$2D$ which does not lie in $S$, which exists by hypothesis. It suffices to
show that we can extend our tuple to a tuple $(M_{\ell})_{\ell\in S\cup T}$
with

$\prod_{\ell\in
T\setminus\\{2\\}}\chi_{\ell}(\mathrm{det}(M_{\ell}))={\mathrm{sgn}}(M_{2})\text{,}$
(4.6)

for then this tuple will lie in the image of $\rho_{E,\Pi^{\prime}}$, where
$\Pi^{\prime}=\prod_{\ell\in S\cup T}\ell$, and so will the original one. But
this is clear: since by hypothesis there is a prime $q\in T\setminus S$, we
can choose $M_{q}$ freely so that (4.6) holds, once all the other $M_{\ell}$
have been prescribed.

We now come to the main part of the proof, namely the case where $D\equiv
1\pmod{4}$ and $S$ contains all prime factors of $2D$. In this case we can
identify $\rho_{E,\Pi}$ with

${\prod_{\ell\in
S}}^{\ast}\mathrm{GL}_{2}(\mathbb{F}_{\ell})=\left\\{(M_{\ell})_{\ell\in
S}:\prod_{\ell\mid
D}\chi_{\ell}(\mathrm{det}(M_{\ell}))={\mathrm{sgn}}(M_{2})\right\\}$

and by Propositions 3.2 and 3.3 and Chebotarev’s density theorem, the desired
density is

$\frac{\\#\left\\{(M_{\ell})_{\ell\in S}\in\prod_{\ell\in
S}^{\ast}\mathrm{GL}_{2}(\mathbb{F}_{\ell}):M_{\ell}\nsim\begin{pmatrix}1&\ast\\\
0&1\end{pmatrix}\right\\}}{\\#\prod_{\ell\in
S}^{\ast}\mathrm{GL}_{2}(\mathbb{F}_{\ell})}\text{.}$

We wish to compare this against

$\delta_{S}=\frac{\\#\left\\{(M_{\ell})_{\ell\in S}\in\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell}):M_{\ell}\nsim\begin{pmatrix}1&\ast\\\
0&1\end{pmatrix}\right\\}}{\\#\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell})}\text{.}$

Since $\prod_{\ell\in S}^{\ast}\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ has index
$2$ in $\prod_{\ell\in S}\mathrm{GL}_{2}(\mathbb{F}_{\ell})$, the ratio
between the desired density and $\delta_{S}$ is

$\frac{2\\#\left\\{(M_{\ell})_{\ell\in S}\in\prod_{\ell\in
S}^{\ast}\mathrm{GL}_{2}(\mathbb{F}_{\ell}):M_{\ell}\nsim\begin{pmatrix}1&\ast\\\
0&1\end{pmatrix}\right\\}}{\\#\left\\{(M_{\ell})_{\ell\in S}\in\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell}):M_{\ell}\nsim\begin{pmatrix}1&\ast\\\
0&1\end{pmatrix}\right\\}}$

and it is clear that each prime $\ell\nmid 2D$ contributes with a factor
$\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ to both the numerator and the
denominator above, so these primes can be ignored, and we are left to compute

$\frac{2\\#\left\\{(M_{\ell})_{\ell\mid 2D}\in\prod_{\ell\in
S}^{\ast}\mathrm{GL}_{2}(\mathbb{F}_{\ell}):M_{\ell}\nsim\begin{pmatrix}1&\ast\\\
0&1\end{pmatrix}\right\\}}{\\#\left\\{(M_{\ell})_{\ell\mid
2D}\in\prod_{\ell\in
S}\mathrm{GL}_{2}(\mathbb{F}_{\ell}):M_{\ell}\nsim\begin{pmatrix}1&\ast\\\
0&1\end{pmatrix}\right\\}}\text{.}$ (4.7)

According to the counting employed in the proof of Lemma 3.1, the denominator
of (4.7) is $\prod_{\ell\mid
2D}(\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})-\ell^{2})$. We now observe that we
can rewrite the numerator as

$\sum\left(1+{\mathrm{sgn}}(M_{2})\prod_{\ell\mid
D}\chi_{\ell}(\mathrm{det}(M_{\ell}))\right)$

where the sum ranges over tuples $(M_{\ell})_{\ell\mid 2D}$ such that each
$M_{\ell}$ is not conjugate to a matrix of the form $\begin{pmatrix}1&\ast\\\
0&1\end{pmatrix}$ for each $\ell$. Indeed, the summand on display equals $2$
if the corresponding tuple is in $\prod_{\ell\in
S}^{\ast}\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ and $0$ otherwise. This equals

$\displaystyle\sum 1+\sum{\mathrm{sgn}}(M_{2})\prod_{\ell\mid
D}\chi_{\ell}(\mathrm{det}(M_{\ell}))$ $\displaystyle=\prod_{\ell\mid
2D}(\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})-\ell^{2})+\left(\sum{\mathrm{sgn}}(M_{2})\right)\prod_{\ell\mid
D}\left(\sum\chi_{\ell}(\mathrm{det}(M_{\ell}))\right)\text{,}$

where the condition on the summands is the same as the one specified before.

One computes directly that $\sum{\mathrm{sgn}}(M_{2})=2$. Moreover, if
$\ell\mid D$ then _all_ matrices in $\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ that
are conjugate to a matrix of the form $\begin{pmatrix}1&\ast\\\
0&1\end{pmatrix}$ have determinant $1$ and as such they lie on the kernel of
$\chi_{\ell}\circ\mathrm{det}$. Therefore, among the matrices considered in
our summation, there are
$\frac{1}{2}\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})-\ell^{2}$ matrices on which
$\chi_{\ell}\circ\mathrm{det}$ evaluates to $1$ and
$\frac{1}{2}\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})$ on which it evaluates to
$-1$. Hence,

$\sum\chi_{\ell}(\mathrm{det}(M_{\ell}))=(\frac{1}{2}\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})-\ell^{2})-\frac{1}{2}\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})=-\ell^{2}\text{.}$

We conclude that (4.7) equals

$\displaystyle\frac{\prod_{\ell\mid
2D}(\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})-\ell^{2})+2\prod_{\ell\mid
D}(-\ell^{2})}{\prod_{\ell\mid
2D}(\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})-\ell^{2})}$
$\displaystyle=1+\frac{2\prod_{\ell\mid D}(-\ell^{2})}{\prod_{\ell\mid
2D}(\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})-\ell^{2})}\text{.}$

Using that
$\\#\mathrm{GL}_{2}(\mathbb{F}_{\ell})=(\ell^{2}-1)(\ell^{2}-\ell)$, this
simplifies to

$1+\frac{\prod_{\ell\mid D}(-\ell^{2})}{\prod_{\ell\mid
D}(\ell^{4}-\ell^{3}-2\ell^{2}+\ell)}=1+\prod_{\ell\mid
D}\frac{-\ell}{\ell^{3}-\ell^{2}-2\ell+1}\text{,}$

as desired. ∎

##  5\. Conclusion

The density we would ultimately like to compute requires a version of Lemma
4.4 where $S$ is allowed to be the set of all primes. We record, however, the
natural conjecture that follows from Lemma 4.4 and that is the expected answer
to our main problem (for Serre curves).

######  Conjecture 5.1.

Define

$C=\prod_{\ell\text{
prime}}\left(1-\frac{\ell}{(\ell-1)^{2}(\ell+1)}\right)\text{.}$

Let $E/\mathbb{Q}$ be a Serre curve, and keep the notation from Proposition
4.3. Then the density of primes $p$ for which
$\mathrm{gcd}(\\#E(\mathbb{F}_{p}),p-1)=1$ is

$\begin{cases}C&\text{ if }D\equiv 0\pmod{4}\\\ \left(1+\prod_{\ell\mid
D}\frac{-\ell}{\ell^{3}-\ell^{2}-2\ell+1}\right)C&\text{ if }D\equiv
1\pmod{4}\text{.}\end{cases}$

It would be tempting to try to establish Conjecture 5.1 by some sort of
limiting procedure, using Lemma 4.4 as an input. However, this does not seem
feasible, and establishing Conjecture 5.1 likely requires bringing more
sophisticated methods into play.

##  6\. Acknowledgements

The author was funded through the Engineering and Physical Sciences Research
Council Doctoral Training Partnership at the University of Warwick. The author
would also like to thank Professor Samir Siksek for suggesting this problem,
as well as for countless support during the elaboration of this manuscript.

## References

  * [1] J. H. Silverman, _The Arithmetic of Elliptic Curves, Second Edition_ , Graduate Texts in Mathematics, Springer 2016
  * [2] S. Siksek, _Integral points on punctured Abelian Varieties_ , arXiv:2105.00570
  * [3] J.-P. Serre, _Propriétés galoisiennes des points d’ordre fini des courbes eliptiques_ , Invent. Math. 15 (1972), no. 4, 259-331
  * [4] N. Jones, _Almost all Elliptic Curves are Serre Curves_ , Transcriptions of the American Mathematical Society, Vol. 362 (2010)
  * [5] T. Barnet-Lamb, D. Geraghty, and T. Gee, _The Sato-Tate conjecture for Hilbert modular forms_ , Journal of the American Mathematical Society 24 (2011), 411–469.
  * [6] T. Barnet-Lamb, D. Geraghty, M. Harris, and R. Taylor, _A family of Calabi-Yau varieties and potential automorphy II_ , Publications of the Research Institute for Mathematical Sciences 47 (2011), 29–98.
  * [7] M. Harris, N. Shepherd-Barron, and R. Taylor, _A family of Calabi-Yau varieties and potential automorphy_ , Annals of Mathematics 171 (2010), 779–813.


# Nucleosome sliding can influence the spreading of histone modifications

Shantanu Kadam1*, Tripti Bameta2*, Ranjith Padinhateeri
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Biosciences and Bioengineering, Indian Institute of Technology
Bombay, Mumbai, India
2Department of Pathology, Tata Memorial Centre, Homi Bhabha National
Institute, Mumbai, India

###### Abstract

Nucleosomes are the fundamental building blocks of chromatin that not only
help in the folding of chromatin but also in carrying epigenetic information.
It is known that nucleosome sliding is responsible for dynamically organizing
chromatin structure and the resulting gene regulation. Since sliding can move
two neighboring nucleosomes physically close or away, can it play a role in
the spreading of histone modifications ? We investigate this by simulating a
stochastic model that couples nucleosome dynamics with the kinetics of histone
modifications. We show that the sliding of nucleosomes can affect the
modification pattern as well as the time it takes to modify a given region of
chromatin. Exploring different nucleosome densities and modification kinetic
parameters, we show that nucleosome sliding can be important for creating
histone modification domains. Our model predicts that nucleosome density
coupled with sliding dynamics can create an asymmetric histone modification
profile around regulatory regions. We also compute the probability
distribution of modified nucleosomes and relaxation kinetics of modifications.
Our predictions are comparable with known experimental results.

##  I Introduction

In cells, DNA is folded and wrapped around octamers of histone proteins
forming an array of nucleosomes. Nucleosomes are considered to be the
fundamental repeating unit of chromatin and its positioning is important for
gene regulation. In typical chromatin, two neighboring nucleosomes are
separated by short segments of linker DNA of lengths ranging from 10 to 60 bp
[1, 2, 3, 4]. Recent advancements in experimental and computational methods
have helped us to understand how nucleosomes are organized along DNA[5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21].

Nucleosomes also carry epigenetic information in the form of histone
modifications apart from the folding of chromatin. Specific amino acid
residues of histone proteins carry chemical modifications like methylation and
acetylations as histone marks[1, 2]. At specified locations on each histone
protein, certain enzymes add or remove relevant chemical groups leading to a
pattern of post-translational modifications along chromatin contour[27, 28].
How these marks get organized along the chromatin is crucial for regulating
cellular processes like gene expression, DNA repair, DNA replication etc [22,
23, 24, 25, 26]. Experimentally, one can measure the pattern of histone
modification at a given instant in a population of cells by ChIP-Seq
methods[32, 33]. However, it is a difficult task to measure the modification
dynamics in individual cells in real-time[29, 30, 31]. Moreover, comprehensive
mechanisms that leads to dynamic histone modification patterns are not fully
understood yet.

The phenomena of spreading and subsequent maintenance of histone modifications
have been experimentally studied with great interest. Several studies[37, 38,
39, 40] have investigated the formation of heterochromatin and epigenetic
inheritance. In these studies, the modified nucleosomes recruit enzymes to
similarly modify neighboring unmodified nucleosomes based on a linear stepwise
process. Also there are studies, where researchers have tried to unravel how
histone modifications spread along chromatin fiber from a given initiation
site[35, 36]. All these studies have considerably contributed to the
understanding of modification spreading.

Several theoretical models[4, 22, 34] have been developed to provide insights
into the dynamics of histone modifications. Over the years, Sneppen and
coworkers have developed models that explain different aspects of histone
modification spreading and inheritance [41, 42, 43, 44]. They have proposed
that long-range interactions lead to a bistable paradigm for a certain range
of parameters. Crabtree et al proposed a linear propagation scheme to explain
patterns in H3K9me3 that involved localized peaks and soft borders of
heterochromatic islands [45]. In their model[45], they incorporated
nucleation, propagation, and turnover rates for modifications, which were
necessary to describe H3K9me3 domains. In a separate work, they extended the
model to estimate several dynamic quantities predicting domain sizes for
different values of rates[46]. This standard model was also extended to
incorporate spreading beyond nearest neighbors. There are also other
stochastic models [47, 48] that include recruitment, diffusion, long- range
interactions leading to the formation of modification patterns. There has been
also a Potts-type model by Zhang et al.[53], stressing the local nature of
interactions, and a model introduced by Binder et al.[54] investigating
epigenetic silencing in eukaryotes. Another set of models explicitly account
for 3D looping and investigate the role of looping in the spreading of histone
modifications [49, 50, 51, 52, 55, 56, 57]. In these models, the coupling of
3D configurations with spreadings of histone modifications is investigated.

Nearly all the modelling studies have assumed that nucleosomes are static and
modifications can spread to the nearest neighbor nucleosomes independent of
the distance between the two nearby nucleosomes. It is plausible that the two
neighboring nucleosomes are far away and the spreading may get hindered due to
the large gap between those nucleosomes[60, 61]. One way two neighboring
nucleosomes can regulate the gap — internucleosomal distance, or linker length
— is via sliding of nucleosomes[58, 59]. None of the existing models account
for the role of the sliding of nucleosomes in the context of modification
spreading. In this manuscript, we propose a stochastic model to study the
spreading and maintenance of histone modifications taking into account the
role of nucleosome sliding.

In the model, we have included sliding of nucleosomes (due to remodeling
complexes), modification of nucleosomes (by modification enzymes), and removal
of a modification mark (by de-modifying enzymes) as kinetic events accompanied
by their respective rates. We aim to explore whether nucleosome sliding events
play any kind of role in spreading the modifications in a particular genomic
region, and how sliding could couple with nucleosome density to determine the
modification dynamics.

We have organized this article as follows: first, we have explained the
features associated with our one-dimensional model along with simulation
details. In the results section, first, we have studied the variation of mean
modification spreading times (MMST) for different sliding rates and de-
modification rates. In the next section, we have computed probabilities of
modified nucleosomes at two different nucleosome densities and compared them
with existing experimental results. We have studied the dynamics of modified
nucleosomes by estimating their statistical quantities. In the last section,
the relaxation dynamics of modified nucleosomes are reported and thoroughly
analyzed their behavior. Finally, we have provided conclusions of our work
along with some suggestions for new experiments to test our predictions.

##  II Model

Figure 1: Schematic of the model: Histone modifying enzymes (“writers”) are
recruited with a rate $k_{recruit}$ at the left boundary. The spreading of the
modification (yellow star) occurs with a rate $k_{spread}$, when the distance
between two neighboring nucleosomes (brown blocks) is less than $l_{c}$. The
de-modification of a nucleosome occurs randomly with a rate $k_{r}$ at any
modified nucleosome. The sliding of a nucleosome occurs at a rate $k_{slide}$
with a step size of 10 bp, provided that there is a free linker DNA space
$\geq l_{c}$ to slide. We take $l_{c}=10$bp in this work.

This section describes the model (see Fig. 1) that we use in our simulations.
We have modelled the spreading of modifications in a quantitative way by
simplifying the 3D structure of DNA and focussing on a small section of DNA
taking it as a 1D lattice. Thus, in the model, the DNA is considered as a one-
dimensional lattice (violet bar) of length $L=5000$ bp. The nucleosomes (brown
rectangular blocks) with indices $i=1,2,...,N$ on the DNA are modeled as hard-
core particles with each one occupying $k$ = 147 bp along the lattice. The
hardcore steric interactions among the nucleosomes are modeled by prohibiting
a lattice site from getting occupied simultaneously by more than one
nucleosome.

In the model, we have considered four kinetic events: i) recruitment of
modification enzymes with the rate $k_{recruit}$, ii) transfer of this
modification enzyme to an unmodified nucleosome with rate $k_{spread}$, iii)
fall of modification enzyme with the rate $k_{r}$ and iv) sliding of
nucleosomes along the DNA with the rate $k_{slide}$. In this study, the
binding and dissociation of nucleosomes are ignored. Hence, we have assumed
that the total number of nucleosomes ($N$) is constant.

We consider a situation where modification enzymes are recruited at a specific
location $i=0$ (left end of the lattice) with a rate $k_{recruit}$ (Fig. 1).
This is implemented when the nucleosome ($i=1$) gets closer to the source of
modification enzymes within a distance (gap between two nucleosomes) less than
$l_{c}$ = 10 bp. The recruited modification enzyme further spreads along the
lattice with the rate $k_{spread}$ to a neighboring unmodified nucleosome
provided that inter-nucleosomal distance is less than $l_{c}$ = 10 bp. One of
the hallmarks of histone modification spreading is the positive feedback[41,
49]. In this model a modified nucleosome inducing modification to spatially
close ($\leq$ 10bp) neighbor is essentially the positive feedback. In the
simulations, it is assumed that the rate of recruitment ($k_{recruit}$) is the
same as the rate of spreading ($k_{spread}$). The modified nucleosome can be
randomly demodified by removal of modification enzyme with a rate $k_{r}$. Our
simulations aim to answer the following question, how nucleosome sliding
affects the spreading of histone modifications ? Hence, we have incorporated
random nucleosome sliding to the left or right with a rate $k_{slide}$ per
nucleosome. This rate represents the rate of reaction by ATP-dependent
chromatin remodelers[62] that are responsible for such repositioning of
nucleosomes in cells. It is assumed that, sliding step size is 10 bp such that
the diffusion constant is $\sim k_{slide}(10{\rm bp})^{2}$. In our studies, it
is assumed that beyond $i=1$ and $i=N$ there are boundary elements, which do
not allow the spreading any further.

We have performed kinetic Monte Carlo simulations using Doob-Gillespie
Algorithm[63, 64, 65, 66] using the rates of sliding, modification, and de-
modification events. All these events are independent of each other. In this
way, the simulation was run for some desired time. The nucleosome density is
calculated as, $\rho=\frac{147*N}{L}$; where $N$ is the total number of
nucleosomes and $L$ is the total length of the 1D lattice. In all the
simulations discussed in this work, the modification rate, $k_{spread}$ is
kept fixed at 1 $s^{-1}$. The rates ($k_{slide}$, $k_{r}$) of all other events
are scaled with $k_{spread}$ giving dimensionless quantities for those
respective rates. The time reported in the simulations is taken in the unit
$\tau_{s}=\frac{1}{k_{spread}}$. All the simulations in this study were
performed by taking an average over 2000 independent runs.

##  III Results

###  III.1 Kinetics of modification spreading

####  III.1.1 Mean Modification Spreading Time: Simulations and Mean Field
Theory

We simulated the spreading of histone modification as discussed in the Model
section for various nucleosome densities computing mean modification spreading
time (MMST) as a function of sliding rates. The MMST is defined as the time
required for the first successful modification of the last nucleosome $(i=N)$,
given that the modification spreads from the initiation site $(i=0)$. These
are similar to the mean first passage time calculations in statistical
mechanics[67], which are a measure of modification spreading time.

We observed that, as the sliding rate increases, the mean modification
spreading time decreases and saturates to a constant value(Fig. 2). For a zero
sliding rate, on a sufficiently long DNA, the probability of finding at least
one pair of nucleosome-neighbors with a gap (linker length) greater than 10 bp
is very high. Hence, for zero sliding rate, the modification may not reach the
other end implying the mean modification spreading time can be infinity. In
order to see the effect of sliding, the mean modification spreading time
(MMST) was computed at three different nucleosome densities (90%, 85%, and
80%). Here, the de-modification rate, $k_{r}$, of nucleosome was kept fixed at
0.01.

Figure 2: Quicker sliding reduces the spreading time: The mean modification
spreading time (MMST) with $k_{slide}$ using simulations (dots) and mean-field
theoretical calculations (lines) for different nucleosome densities, viz. :
90% (red), 85% (green) and 80% (blue). All rates are measured in units of
$k_{spread}$.

We also used a mean-field theoretical study to understand how MMST would vary
with sliding rates for different densities. For this calculation, at a given
density, we have taken nucleosomes to be homogeneously distributed along the
lattice. The effect of nucleosome sliding was incorporated in the effective
spreading rate of modification, $k_{se}$, which is a function of the sliding
rate of the nucleosome and average inter-nucleosomal distance. The mean
modification spreading time ($T_{i\rightarrow n}$), from $i^{th}$ to $n^{th}$
nucleosome follows the difference equation[68]:

$T_{i\to n}=\dfrac{1}{k_{\rm se}+k_{r}}+\dfrac{k_{se}}{k_{se}+k_{r}}T_{i+1\to
n}+\dfrac{k_{r}}{k_{se}+k_{r}}T_{i-1\to n}$ (1)

where index $i$ varies from $0$ to $N$. The index $i=0$ represent nucleation
site, while ($1$) and ($N$) are indices of first and last nucleosome
respectively. The $N$ is the total number of nucleosomes in the model. From
solving this $N+1$ set of linear equations (see Appendix) :

$T_{0\to N}=\frac{1}{{(k_{\rm se})^{N}}}\sum_{\ell=1}^{N}(N-\ell+1)(k_{\rm
se})^{N-\ell}~{}~{}k_{\rm r}^{\ell-1}$ (2)

It is expected, when de-modification rate of nucleosomes is set to zero (i.e.
$k_{\rm r}=0$) (2) reduces to $T_{0\to N}=\frac{N}{k_{\rm se}}$. In this
equation, $k_{\rm se}$ depends on the sliding rate of nucleosome and the
inter-nucleosomal distance via the relation[69, 70]:

$k_{se}=\frac{l_{s}^{2}\times k_{slide}}{gap^{2}}$ (3)

where $l_{s}=10$ bp is the step size of sliding events, and the gap is the
linker length. This relation can be understood as an inverse of the time scale
of the meeting of two nucleosomes. In Fig. 2, we plot the mean-field theory
results (Eq.2, curves) along with simulation results (dots). For certain
nucleosome densities, both results are comparable. A significant variation was
observed in MMST when density was varied from 90% to 85%. The higher values of
MMST at lower densities were a signature of the presence of long gaps (greater
than 10 bp) between nucleosomes, which eventually slowdowns the spreading of
modifications. This also implies, when $k_{slide}\ll k_{spread}$ MMST is
higher, while in the opposite limit it is found to be smaller.

We have got elevated profiles of MMST at 85% and 80% nucleosome densities
implying reduced densities contributing towards an increase in MMST values.
For low sliding rates (less than or around 0.5), the MMST changes a lot with
the sliding rate indicating that the longer times it takes to spread the
modifications across the lattice. For higher sliding rates, MMST change is
small. Altogether, all these results imply that the sliding of nucleosomes can
play a significant role in spreading the modification across the lattice.

####  III.1.2 De-modification Events Increase Spreading Time

In reality, modified nucleosomes can get de-modified and such de-modification
events may be crucial in some contexts. To examine the effect of de-
modification on MMST, we varied the de-modification rate ($k_{r}$) over a
range. In Fig. 3 we present our results for MMST by changing the de-
modification rates of nucleosomes for four different sliding rates 0.5 (red),
1 (green), 2 (blue) and 6 (magenta). Since, the residence time of the
modification[47, 71] is more than spreading time, we take de-modification rate
$k_{r}\ll 1$ in units of $\frac{1}{\tau_{s}}$.

Figure 3: Spreading time increases with de-modification events: The mean
modification spreading time (MMST) as a function of de-modification rate
$k_{r}$ for different nucleosome densities : (A) 90%, (B) 85% and (C) 80%.
Different curves are for sliding rates 0.5, 1, 2, and 6 from top (red) to
bottom (magenta). All rates are measured in units of $k_{spread}$.

In Figs. 3 (A), (B) and (C) MMST results are plotted for different nucleosome
densities 90%, 85%, and 80% respectively. Here, one sees an interplay between
nucleosome de-modification rate and sliding rate. An increase in de-
modification rates contributes to the corresponding increase in the MMST. It
would take a longer time for spreading the modification across the lattice for
larger values of de-modification rates. However, an increase in sliding rates
from 0.5 to 6 (red, green, blue, magenta curves) has had contributed to a
substantial decrease in the MMST. It was found that with a reduction in
nucleosome density, the gaps between nucleosomes as well as de-modification
rates have contributed to an increase in modification spreading times.

###  III.2 Domains of modified nucleosomes : effect of sliding and de-
modification events

In this section, we discuss how the modified domains of nucleosomes are
maintained for different sliding and de-modification rates at a fixed
nucleosome density. We have calculated the probability of modified nucleosomes
at a steady state. The simulations were carried out by varying sliding rates
whilst fixing $k_{r}$ = 0.1.

Figure 4: Probabilities with sliding and de-modification rates: The
probability of modified nucleosomes for sliding and de-modification rates at
80% nucleosome density: (A) for $k_{r}$ = 0.1, $k_{slide}$ = 2.2 (red), 1
(magenta), 0.5 (green) and 0.1 (blue) and (B) for $k_{slide}$ = 0.25, $k_{r}$
= 0.02 (red), 0.03 (green), 0.05 (blue), 0.1 (magenta). The black squares
depict experimental data (Hathaway et al. ([45])) for MEF cells and grey
squares for ES cells.(C) Probabilities at 80% density (left of initiation
site) for $k_{slide}$ = 0.1 (red), 1 (green) and 2.2 (blue) open circles and
90% density (right of initiation site) for $k_{slide}$ = 0.1 (red) and 0.075
(blue) filled circles. All times are measured in units of $\tau_{s}$.

The probability of finding modified nucleosomes at any location along the DNA
contour, at 80% nucleosome density, is depicted in Fig. 4 with different
colors representing different sliding and de-modification rates. The
probability sharply decays as we decrease the sliding rate (Fig. 4(A)). In
other words, the modification pattern tends to get more localized with a peak
at the initiation (source) site. This sharp decay in modification profile is
similar to what is observed in experiments near the nucleation site. For
$k_{slide}$ = 1, our simulation results are comparable to experimental H3K9me3
ChIP data from Hathaway et al.[45] of mouse embryonic fibroblasts (MEFs). A
similar set of simulations were repeated at 90% nucleosome density for
different sliding rates (see Fig S4-A). In this case, density has played an
important role in modification of these nucleosomes along with their sliding
rates. A significant decrease in probabilities at the end of lattice is a
signature of de-modification events dominating over the spreading of
modifications.

As a next step, simulations were carried out by varying de-modification rates
while fixing the $k_{slide}$ = 0.25 (Fig. 4(B)). For small $k_{r}$ values,
like 0.02 and 0.03, the modification pattern has a larger spread. It was found
that for $k_{r}$ = 0.02, our simulation results are comparable to experimental
H3K9me3 ChIP data from Hathaway et al.[45] of Embryonic Stem (ES) cells. For
relatively higher $k_{r}$ values localized modification pattern is observed
with a peak near the initiation site. At $k_{r}$ = 0.05 the simulated
modification profile is comparable to what is observed in experiments for MEF
cells.

Our results so far suggest that combination of different sliding rates and
nucleosome density can lead to very different modification patterns. Since it
is plausible that on either side of certain boundary or boundary elements
(e.g. transcription start site), nucleosome density and action of chromatin
remodellers could be very different. We explore this to examine whether this
can lead to an asymmetry in nucleosome modification spread patterns[55]. We
simulated nucleosomes with different densities and sliding rates on either
side of a boundary (initiation site). The results are in Fig. 4(C), where the
asymmetry in the positioning of modified nucleosomes about the initiation site
can be seen. The density of the nucleosomes are always fixed at 90% on the
right-hand side of the initiation site and at 80% on the left-hand side. Then
different combinations of sliding rates are taken on both sides of initiation
site. The top blue curve (with open circles on the left of initiation site for
$k_{slide}$ = 2.2 and filled circles on right of initiation site for
$k_{slide}$ = 0.075) has least asymmetry; the bottom red curve (with open
circles on the left of initiation site and filled circles on right of
initiation site for $k_{slide}$ = 0.1) has maximum asymmetry. Thus, we suggest
that variability in sliding rate and nucleosome densities is a potential way
of explaining the asymmetry across various boundary elements.

These results together show that the interplay between nucleosome de-
modification rate and the sliding rate determines the profile of the
modification pattern. At low sliding and high de-modification rates, the
modification is peaked near the source and decays quickly. For high sliding
rate and low de-modification rate, the decay is more gradual.

###  III.3 Dynamics of modified nucleosomes : estimation of statistical
quantities

In this subsection, we discuss the time evolutions in the amount of
modification and its fluctuations at biologically relevant nucleosome
densities (80% and 90%). For each time step, we plot the mean fraction of
modified nucleosomes $\frac{\rm{N_{\rm{m}}}}{\rm{N}}$, where $\rm{N_{\rm{m}}}$
is the mean number of modified nucleosomes (Fig. 5). The fluctuations were
quantified using the coefficient of variation (CV), a dimensionless quantity
defined as the ratio of standard deviation to mean. These quantities were
calculated for different sliding rates $k_{slide}$ = 0.5, 1, 2. $k_{r}$ was
kept fixed at 0.1. At $t=0$, $\rm{N_{\rm{m}}}$ is taken as zero. As expected,
$\frac{\rm{N_{\rm{m}}}}{\rm{N}}$ increases with time and saturates. Here
different trajectories (colors) correspond to different rates of sliding. At
higher density, the variation among $\frac{\rm{N_{\rm{m}}}}{\rm{N}}$ for
different sliding rates (i.e. different curves with different colors) is
small. However, at lower density (80%), the mean dynamics show huge variation
as we change the sliding rates, suggesting that steady-state as well as the
dynamics at low densities are crucially affected by the rates.

Figure 5: Statistics of modified nucleosomes: The (A) fraction of modified
nucleosomes and (B) coefficient of variation (CV) of modified nucleosomes for
90% nucleosome density (N = 31 nucleosomes) (C) fraction of modified
nucleosomes and (D) coefficient of variation (CV) of modified nucleosomes for
80% nucleosome density (N = 27 nucleosomes). All plots are for different
nucleosome sliding rates $k_{slide}$ = 0.5 (red), 1 (green) and 2 (blue) for
lattice length $L$ = 5000 bp. All times are measured in units of $\tau_{s}$.
Figure 6: Number distribution profiles of modified and de-modified
nucleosomes: The probability distributions of number of (A) modified and (B)
de-modified nucleosomes at 90% nucleosome density. The probability
distributions of number of (C) modified and (D) de-modified nucleosomes at 80%
nucleosome density for nucleosome sliding rates $k_{slide}$ = 0.5 (red), 1
(green) and 2 (blue). All rates are measured in units of $k_{spread}$.

We find that CV values are decreasing as a function of time. This could be
because at early times the modification numbers are less and hence the
fluctuation is high. At 90% density, CV is smaller than unity (standard
deviation smaller than the mean) and CV is similar at the steady-state for all
sliding rates. However, for 80% nucleosome density, the CVs show well-
separated steady states, and fluctuations are comparable or bigger than the
mean (CV is comparable or above $1$). We found almost the same steady states
for fraction of modified nucleosomes and CV, when the simulation was done at
80% nucleosome density with all nucleosomes modified (N = 27) at t = 0. It
implies that the final states are independent of the initial conditions.
(whether all nucleosomes or none of them were modified) (see Fig S1)

The probability distributions for the number of modified and de-modified
nucleosomes at 90% and 80% nucleosome densities are depicted in Fig. 6. These
distributions were obtained by finding nucleosome numbers at steady states. In
Fig. 6(A), for given sliding rates distributions of modified nucleosomes are
negatively skewed. It was found that at 90% density almost all nucleosomes
have got modified for lower (0.5) and higher (2) values of sliding rates. This
implies at high nucleosome density modification is not much affected by
sliding rates of nucleosomes. The distribution of de-modified nucleosomes is
shown in Fig. 6(B). At 80% nucleosome density, in Fig. 6(C and D)
qualitatively same distributions were observed for low and high values of
sliding rates. At a lower sliding rate (0.5) less number of nucleosomes were
modified; while at a higher sliding rate modification spreads along the
lattice modifying a large number of nucleosomes.

###  III.4 Relaxation dynamics of modifications when the initiation site is
removed

Figure 7: Relaxation of modified nucleosomes for different sliding rates : The
time profiles of fraction of modified nucleosomes at nucleosome densities (A)
90% (N = 31 nucleosomes) and (B) 85% (N = 29 nucleosomes) when all nucleosomes
are modified at t = 0 and $k_{recruit}=0$. At nucleosome densities (C) 90% and
(D) 85% at t = 0, $k_{recruit}=1$ and all nucleosomes are unmodified. At
steady state (black vertical line) of modification we take $k_{recruit}=0$.
All the results are for sliding rates $k_{slide}$ = 0.1 (red), 0.2 (green),
0.5 (blue). The fitted curve is shown in black and all times are measured in
units of $\tau_{s}$.

The cell actively maintains the average number of modified nucleosomes by
constantly inserting modifications. In this subsection, we will examine how
the number of modified nucleosomes decreases as a function of time if the
nucleation site removed (spreading from the nucleation site is switched off).
We will discuss two different cases. In the first case, at t = 0 all
nucleosomes are assumed to be in the modified state. As the time starts, we
stop influx of histone modification enzymes from the initiation site
($k_{recruit}=0$). That is, starting from a fully modified state, we fix
$k_{recruit}=0$ and simulate dynamics taking all other events. In the second
case, we simulate the full system to obtain the steady state, and at some time
point in steady state, we set $k_{recruit}=0$.

In Fig. 7(A) we present the results for the first case with higher density
(90%) and different nucleosome sliding rates (0.2 (green), 0.5 (blue) and (0.1
(red)). We observe stable patterns of modified nucleosomes are maintained for
sliding rates 0.2 and 0.5 in respective steady states. However, when the
sliding rate is low (0.1 (red)) the fraction of modified nucleosomes decreases
as a function of time. This curve was fitted with an exponentially decaying
function giving rise to a decay rate of the order of $10^{-6}$ per unit time.
This is analogous to the effective de-modification rate in our simulations.
The loss of steady-state pattern emphasizes the importance of sliding events
even at such a higher density. However, at lower densities in (85%) Fig. 7(B),
such steady-state patterns are absent for all the sliding rates simulated.
However, for higher sliding rate, the mean modification decays slowly.

In the second case, we simulated the full system (taking all events discussed
in subsection C) until the steady state. At a particular time point indicated
by the vertical bar in Fig. 7(C) and (D), influx of modification was switched
off ($k_{recruit}$ is set to $0$). For higher sliding rates (0.2 (green) and
0.5 (blue)) and higher density, there is no dip in the fraction of modified
nucleosomes. However, for a lower sliding rate (0.1(red)) the fraction of
modified nucleosomes decreases a bit but stays near the steady state. At 85%
nucleosome density, the fraction of modified nucleosomes decay by slightly but
stays closer to their respective steady states. Fitting the decay curve with
an exponential function, the decay constant at 85% nucleosome density is found
to be around $10^{-3}$ per unit time (see Supplementary Table S1 and Table
S2).

In short, the interplay between density and sliding decide the relaxation
dynamics of modifications; an increase in $k_{slide}$ values (from 0.1 to 0.5)
the modified nucleosomes decay slowly. The sliding of the nucleosome is
contributing to slowing down the decay (also see Supporting information Fig S2
and Fig S3).

##  IV Conclusions and suggestions for experiments

In this article, we have discussed how nucleosome sliding may affect the
spreading of histone modifications along the lattice. The spread of
modifications is quantified by estimating the mean modification spreading time
(MMST). It was found that for low nucleosomal sliding rates it takes a longer
time for spreading the modification across the lattice; while it took less
time for higher sliding rates. We have confirmed these findings by doing an
analytical estimation of MMST using a mean field theory. It was also found
that the interplay between nucleosome sliding events and nucleosome density
determines spreading times. The larger de-modification rates have contributed
to enhanced modification spreading times, but the sliding rates of nucleosomes
have helped to restrict them. The dynamics of modified nucleosomes were
studied by computing statistical quantities like fluctuations and probability
distributions, which were found to be dependent upon nucleosome densities. We
show that for certain densities and sliding rates, the nucleosome modification
pattern is localised in a region, as seen in experiments. This work also shows
that certain parameters can give raise to asymmetric nature of nucleosome
modifications about the initiation site. The interplay between sliding events
and density of nucleosomes also influence the relaxation dynamics of
modifications. Overall, the proposed model gives insights into the role of
sliding events and how the interplay between density and sliding can be an
important determinant.

The main prediction of the paper is how the sliding of nucleosomes and
nucleosome densities can influence the spreading of modifications. This can be
tested by developing appropriate mutants of nucleosome sliding enzymes and
examining whether the mutations affect the spreading of modifications or not.
One may also design different chromatin arrays having very different
nucleosome densities (nucleosome repeat lengths) and examine how these would
influence the spreading of modifications.

## Acknowledgments

This research has been supported by funds to SK from Department of Science and
Technology, India under Science and Engineering Research Board, National Post-
doctoral Fellowship (NPDF) with file number: PDF/2017/002502 and to RP from
Department of Biotechnology, India Grant BT/HRD/NBA/39/12/2018-19.

##  Appendix A Recursive relation of mean first passage time from survival
probability

A simplified mean field calculation of this model can give a closed form
expression for MMST as a function of rates. For this calculation, we take
nucleosomes homogeneously distributed along the lattice. The modification
spreads from the nucleation site with an effective rate $k_{se}$, which
depends on sliding rate and inter-nucleosomal distance by the following
relation[69, 70]:

$\displaystyle k_{se}=\frac{l_{s}^{2}\times k_{slide}}{gap^{2}}$ (4)

For MMST calculations, we first calculate survival probability ($S_{i\to
n}(t)$) which is defined as: the probability that modification has not reached
to the $n^{th}$ nucleosome till time $t$ given that at $t=0$ $i^{th}$
nucleosome was already modified[68].

$\displaystyle\dfrac{\partial S_{i\to n}(t)}{\partial t}$ $\displaystyle=$
$\displaystyle k_{se}S_{i+1\to n}(t)+k_{r}S_{i-1\to
n}-\left(k_{se}+k_{r}\right)S_{i\to n}(t).$ (5)

We can write first passage time $T_{i\to n}$ as;

$\displaystyle T_{i\to n}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}t~{}F_{i\to n}(t)dt$ (6) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}t\dfrac{-\partial S_{i\to n}(t)}{\partial t}dt$
$\displaystyle=$ $\displaystyle-t~{}S_{i\to
n}(t)|_{0}^{\infty}-\int_{0}^{\infty}1.\left(-S_{i\to n}(t)\right)dt$
$\displaystyle=$ $\displaystyle\int_{0}^{\infty}S_{i\to n}(t)dt$ (7)

First term in Eq.(A3) vanishes because for any bounded survival probability at
large time is $\approx 0$, and it approaches zero much faster than $1/t$.

$\displaystyle\int_{0}^{\infty}\dfrac{\partial S_{i\to n}(t)}{\partial t}dt$
$\displaystyle=$ $\displaystyle k_{se}\int_{0}^{\infty}S_{i+1\to
n}(t)dt+k_{r}\int_{0}^{\infty}S_{i-1\to
n}dt-\left(k_{se}+k_{r}\right)\int_{0}^{\infty}S_{i\to n}(t)dt$ $\displaystyle
S_{i\to n}(\infty)-S_{i\to n}(0)$ $\displaystyle=$ $\displaystyle
k_{se}T_{i+1\to n}+k_{r}T_{i-1\to n}-\left(k_{se}+k_{r}\right)T_{i\to n}$ (8)
$\displaystyle T_{i\to n}$ $\displaystyle=$
$\displaystyle\dfrac{1}{k_{se}+k_{r}}+\dfrac{k_{r}}{k_{se}+k_{r}}T_{i-1\to
n}+\dfrac{k_{se}}{k_{se}+k_{r}}T_{i+1\to n}$ (9)

In Eq.(A5), we have used two properties of survival probability that a system
should survive with probability $1$ at $t=0$ and modification should survive
with probability $\approx 0$ at very large time.

Now using Eq.(A6), we can write a recursive equation between MMSTs : ($T_{i\to
N}\rightarrow{\rm MMST~{}from}~{}{i}^{\rm th}~{}{\rm
nucleosome~{}to~{}}{N}^{\rm th}~{}{\rm nucleosome}$)

$\displaystyle T_{0\to N}$ $\displaystyle=$
$\displaystyle\dfrac{1}{k_{recruit}}+T_{1\to N}$ $\displaystyle T_{1\to N}$
$\displaystyle=$
$\displaystyle\dfrac{1}{k_{se}+k_{r}}+\dfrac{k_{se}}{k_{se}+k_{r}}T_{2\to
N}+\dfrac{k_{r}}{k_{se}+k_{r}}T_{0\to N}$ $\displaystyle T_{2\to N}$
$\displaystyle=$
$\displaystyle\dfrac{1}{k_{se}+k_{r}}+\dfrac{k_{se}}{k_{se}+k_{r}}T_{3\to
N}+\dfrac{k_{r}}{k_{se}+k_{r}}T_{1\to N}$ . . . $\displaystyle T_{i\to N}$
$\displaystyle=$
$\displaystyle\dfrac{1}{k_{se}+k_{r}}+\dfrac{k_{se}}{k_{se}+k_{r}}T_{i+1\to
N}+\dfrac{k_{r}}{k_{se}+k_{r}}T_{i-1\to N}$ . . . $\displaystyle T_{N-1\to N}$
$\displaystyle=$
$\displaystyle\dfrac{1}{k_{se}+k_{r}}+\dfrac{k_{r}}{k_{se}+k_{r}}T_{N-2\to N}$
$\displaystyle T_{N\to N}$ $\displaystyle=$ $\displaystyle 0$ (10)

We can write these equations in matrix form:

$\displaystyle\mathbb{M}1~{}\boldsymbol{T}=\boldsymbol{B}$
$\displaystyle\begin{bmatrix}k_{\rm recruit}&-k_{\rm recruit}&&&0\\\ -k_{\rm
r}&k_{\rm se}+k_{\rm r}&-k_{\rm se}&\ddots&\ddots\\\
0&\ddots&\ddots&\ddots&\ddots\\\ &&&&\\\ &\ddots&-k_{\rm r}&k_{\rm se}+k_{\rm
r}&-k_{\rm se}\\\ 0&&&-k_{\rm r}&k_{\rm se}+k_{\rm
r}\end{bmatrix}\begin{bmatrix}T_{0\to N}\\\ T_{1\to N}\\\ .\\\ .\\\ .\\\ .\\\
T_{N-2\to N}\\\ T_{N-1\to N}\end{bmatrix}=\begin{bmatrix}1\\\ 1\\\ .\\\ .\\\
.\\\ .\\\ .\\\ 1\end{bmatrix}$ (11)

$k_{\rm recruit}$ and $k_{\rm se}$ are function of nucleosome density and
$k_{\rm slide}$ We solve the above equation by taking inverse of matrix
$\mathbb{M}$ to get $T_{1\to N}$. For simplicity first take $k_{\rm
recruit}=k_{\rm se}$

$\displaystyle T_{0\to N}=\dfrac{1}{{(k_{\rm
se})^{N}}}\sum_{\ell=1}^{N}(N-\ell+1)k_{\rm se}^{N-\ell}~{}~{}k_{\rm
r}^{\ell-1}$ (12)

## References

  * [1] Alberts B.  Molecular Biology of the Cell: Hauptbd.  Garland; 2002.
  * [2] Van Holde KE.  Chromatin.  Springer Science & Business Media; 2012.
  * [3] Kornberg RD.  Chromatin structure: a repeating unit of histones and DNA.  Science. 1974;184(4139):868–871.
  * [4] Cortini R, Barbi M, Caré BR, Lavelle C, Lesne A, Mozziconacci J, et al.  The physics of epigenetics.  Reviews of Modern Physics. 2016;88(2):025002.
  * [5] Kaplan N, Moore IK, Fondufe-Mittendorf Y, Gossett AJ, Tillo D, Field Y, et al.  The DNA-encoded nucleosome organization of a eukaryotic genome.  Nature. 2009;458(7236):362–366.
  * [6] Brogaard K, Xi L, Wang JP, Widom J.  A map of nucleosome positions in yeast at base-pair resolution.  Nature. 2012;486(7404):496–501.
  * [7] Raveh-Sadka T, Levo M, Shabi U, Shany B, Keren L, Lotan-Pompan M, et al.  Manipulating nucleosome disfavoring sequences allows fine-tune regulation of gene expression in yeast.  Nature Genetics. 2012;44(7):743.
  * [8] Milani P, Chevereau G, Vaillant C, Audit B, Haftek-Terreau Z, Marilley M, et al.  Nucleosome positioning by genomic excluding-energy barriers.  Proceedings of the National Academy of Sciences. 2009;106(52):22257–22262.
  * [9] Zhang Z, Wippo CJ, Wal M, Ward E, Korber P, Pugh BF.  A packing mechanism for nucleosome organization reconstituted across a eukaryotic genome.  Science. 2011;332(6032):977–980.
  * [10] Kornberg RD, Stryer L.  Statistical distributions of nucleosomes: nonrandom locations by a stochastic mechanism.  Nucleic Acids Research. 1988;16(14):6677–6690.
  * [11] Padinhateeri R, Marko JF.  Nucleosome positioning in a model of active chromatin remodeling enzymes.  Proceedings of the National Academy of Sciences. 2011;108(19):7799–7803.
  * [12] Parmar JJ, Marko JF, Padinhateeri R.  Nucleosome positioning and kinetics near transcription-start-site barriers are controlled by interplay between active remodeling and DNA sequence.  Nucleic Acids Research. 2014;42(1):128–136.
  * [13] Teif VB, Rippe K.  Predicting nucleosome positions on the DNA: combining intrinsic sequence preferences and remodeler activities.  Nucleic Acids Research. 2009;37(17):5641–5655.
  * [14] van der Heijden T, van Vugt JJ, Logie C, van Noort J.  Sequence-based prediction of single nucleosome positioning and genome-wide nucleosome occupancy.  Proceedings of the National Academy of Sciences. 2012;109(38):E2514–E2522.
  * [15] Morozov AV, Fortney K, Gaykalova DA, Studitsky VM, Widom J, Siggia ED.  Using DNA mechanics to predict in vitro nucleosome positions and formation energies.  Nucleic Acids Research. 2009;37(14):4707–4722.
  * [16] Chereji RV, Ramachandran S, Bryson TD, Henikoff S.  Precise genome-wide mapping of single nucleosomes and linkers in vivo.  Genome Biology. 2018;19(1):1–20.
  * [17] Jiang Z, Zhang B.  Theory of active chromatin remodeling.  Physical Review Letters. 2019;123(20):208102.
  * [18] Jiang Z, Zhang B.  On the role of transcription in positioning nucleosomes.  PLoS Computational Biology. 2021;17(1):e1008556.
  * [19] Teif VB.  Nucleosome positioning: resources and tools online.  Briefings in bioinformatics. 2016;17(5):745–757.
  * [20] Teif VB, Clarkson CT. Nucleosome Positioning.; 2019.
  * [21] Bameta T, Das D, Padinhateeri R.  Coupling of replisome movement with nucleosome dynamics can contribute to the parent–daughter information transfer.  Nucleic Acids Research. 2018;46(10):4991–5000.
  * [22] Erdel F.  How communication between nucleosomes enables spreading and epigenetic memory of histone modifications.  BioEssays. 2017;39(12):1700053.
  * [23] Allis CD, Jenuwein T, Reinberg D.  Overview and concepts.  Epigenetics. 2007;1.
  * [24] Ramachandran S, Henikoff S.  Replicating nucleosomes.  Science Advances. 2015;1(7):e1500587.
  * [25] Margueron R, Justin N, Ohno K, Sharpe ML, Son J, Drury Iii WJ, et al.  Role of the polycomb protein EED in the propagation of repressive histone marks.  Nature. 2009;461(7265):762–767.
  * [26] Bonasio R, Tu S, Reinberg D.  Molecular signals of epigenetic states.  Science. 2010;330(6004):612–616.
  * [27] Pokholok DK, Harbison CT, Levine S, Cole M, Hannett NM, Lee TI, et al.  Genome-wide map of nucleosome acetylation and methylation in yeast.  Cell. 2005;122(4):517–527.
  * [28] Lien WH, Guo X, Polak L, Lawton LN, Young RA, Zheng D, et al.  Genome-wide maps of histone modifications unwind in vivo chromatin states of the hair follicle lineage.  Cell Stem Cell. 2011;9(3):219–232.
  * [29] Hayashi-Takanaka Y, Yamagata K, Nozaki N, Kimura H.  Visualizing histone modifications in living cells: spatiotemporal dynamics of H3 phosphorylation during interphase.  Journal of Cell Biology. 2009;187(6):781–790.
  * [30] Hayashi-Takanaka Y, Yamagata K, Wakayama T, Stasevich TJ, Kainuma T, Tsurimoto T, et al.  Tracking epigenetic histone modifications in single cells using Fab-based live endogenous modification labeling.  Nucleic Acids Research. 2011;39(15):6475–6488.
  * [31] Sato Y, Hilbert L, Oda H, Wan Y, Heddleston JM, Chew TL, et al.  Histone H3K27 acetylation precedes active transcription during zebrafish zygotic genome activation as revealed by live-cell analysis.  Development. 2019;146(19).
  * [32] Zhang B, Wu X, Zhang W, Shen W, Sun Q, Liu K, et al.  Widespread enhancer dememorization and promoter priming during parental-to-zygotic transition.  Molecular Cell. 2018;72(4):673–686.
  * [33] Lindeman LC, Winata CL, Aanes H, Mathavan S, Aleström P, Collas P.  Chromatin states of developmentally-regulated genes revealed by DNA and histone methylation patterns in zebrafish embryos.  International Journal of Developmental Biology. 2010;54(5):803–813.
  * [34] Zhang T, Cooper S, Brockdorff N.  The interplay of histone modifications–writers that read.  EMBO Reports. 2015;16(11):1467–1481.
  * [35] Allshire RC, Madhani HD.  Ten principles of heterochromatin formation and function.  Nature Reviews Molecular Cell Biology. 2018;19(4):229–244.
  * [36] Buscaino A, Lejeune E, Audergon P, Hamilton G, Pidoux A, Allshire RC.  Distinct roles for Sir2 and RNAi in centromeric heterochromatin nucleation, spreading and maintenance.  The EMBO Journal. 2013;32(9):1250–1264.
  * [37] Hall IM, Shankaranarayana GD, Noma Ki, Ayoub N, Cohen A, Grewal SI.  Establishment and maintenance of a heterochromatin domain.  Science. 2002;297(5590):2232–2237.
  * [38] Schotta G, Ebert A, Krauss V, Fischer A, Hoffmann J, Rea S, et al.  Central role of Drosophila SU (VAR) 3–9 in histone H3-K9 methylation and heterochromatic gene silencing.  The EMBO Journal. 2002;21(5):1121–1131.
  * [39] Felsenfeld G, Groudine M.  Controlling the double helix.  Nature. 2003;421(6921):448–453.
  * [40] Zhang K, Mosch K, Fischle W, Grewal SI.  Roles of the Clr4 methyltransferase complex in nucleation, spreading and maintenance of heterochromatin.  Nature Structural & Molecular Biology. 2008;15(4):381–388.
  * [41] Dodd IB, Micheelsen MA, Sneppen K, Thon G.  Theoretical analysis of epigenetic cell memory by nucleosome modification.  Cell. 2007;129(4):813–822.
  * [42] Sneppen K, Micheelsen MA, Dodd IB.  Ultrasensitive gene regulation by positive feedback loops in nucleosome modification.  Molecular Systems Biology. 2008;4(1):182.
  * [43] Dodd IB, Sneppen K.  Barriers and silencers: a theoretical toolkit for control and containment of nucleosome-based epigenetic states.  Journal of Molecular Biology. 2011;414(4):624–637.
  * [44] Obersriebnig MJ, Pallesen EM, Sneppen K, Trusina A, Thon G.  Nucleation and spreading of a heterochromatic domain in fission yeast.  Nature Communications. 2016;7(1):1–11.
  * [45] Hathaway NA, Bell O, Hodges C, Miller EL, Neel DS, Crabtree GR.  Dynamics and memory of heterochromatin in living cells.  Cell. 2012;149(7):1447–1460.
  * [46] Hodges C, Crabtree GR.  Dynamics of inherently bounded histone modification domains.  Proceedings of the National Academy of Sciences. 2012;109(33):13296–13301.
  * [47] Anink-Groenen LC, Maarleveld TR, Verschure PJ, Bruggeman FJ.  Mechanistic stochastic model of histone modification pattern formation.  Epigenetics & chromatin. 2014;7(1):30.
  * [48] Ancona M, Michieletto D, Marenduzzo D.  Competition between local erasure and long-range spreading of a single biochemical mark leads to epigenetic bistability.  Phys Rev E. 2020;101:042408.  doi:10.1103/PhysRevE.101.042408.
  * [49] Erdel F, Greene EC.  Generalized nucleation and looping model for epigenetic memory of histone modifications.  Proceedings of the National Academy of Sciences. 2016;113(29):E4180–E4189.
  * [50] Michieletto D, Orlandini E, Marenduzzo D.  Polymer model with epigenetic recoloring reveals a pathway for the de novo establishment and 3D organization of chromatin domains.  Physical Review X. 2016;6(4):041047.
  * [51] Jost D, Vaillant C.  Epigenomics in 3D: importance of long-range spreading and specific interactions in epigenomic maintenance.  Nucleic Acids Research. 2018;46(5):2252–2264.
  * [52] Sandholtz SH, Beltran BG, Spakowitz AJ.  Physical modeling of the spreading of epigenetic modifications through transient DNA looping.  Journal of Physics A: Mathematical and Theoretical. 2019;52(43):434001.
  * [53] Zhang H, Tian XJ, Mukhopadhyay A, Kim K, Xing J.  Statistical mechanics model for the dynamics of collective epigenetic histone modification.  Physical review letters. 2014;112(6):068101.
  * [54] Binder H, Steiner L, Rohlf T, Prohaska S, Galle J.  Transcriptional memory emerges from cooperative histone modifications.  Nature Precedings. 2011; p. 1–1.
  * [55] Sandholtz SH, Kannan D, Beltran BG, Spakowitz AJ.  Chromosome Structural Mechanics Dictates the Local Spreading of Epigenetic Marks.  Biophysical Journal. 2020;119(8):1630–1639.
  * [56] Sandholtz SH, MacPherson Q, Spakowitz AJ.  Physical modeling of the heritability and maintenance of epigenetic modifications.  Proceedings of the National Academy of Sciences. 2020;117(34):20423–20429.
  * [57] Katava M, Shi G, Thirumalai D.  Chromatin dynamics controls epigenetic domain formation.  bioRxiv. 2021;.
  * [58] Narlikar GJ, Sundaramoorthy R, Owen-Hughes T.  Mechanisms and functions of ATP-dependent chromatin-remodeling enzymes.  Cell. 2013;154(3):490–503.
  * [59] Zhou CY, Johnson SL, Gamarra NI, Narlikar GJ.  Mechanisms of ATP-dependent chromatin remodeling motors.  Annual Review of Biophysics. 2016;45:153–181.
  * [60] Szerlong HJ, Hansen JC.  Nucleosome distribution and linker DNA: connecting nuclear function to dynamic chromatin structure.  Biochemistry and Cell Biology. 2011;89(1):24–34.
  * [61] Bednar J, Horowitz RA, Grigoryev SA, Carruthers LM, Hansen JC, Koster AJ, et al.  Nucleosomes, linker DNA, and linker histone form a unique structural motif that directs the higher-order folding and compaction of chromatin.  Proceedings of the National Academy of Sciences. 1998;95(24):14173–14178.
  * [62] Swygert SG, Peterson CL.  Chromatin dynamics: interplay between remodeling enzymes and histone modifications.  Biochimica et Biophysica Acta (BBA)-Gene Regulatory Mechanisms. 2014;1839(8):728–736.
  * [63] Doob JL.  Topics in the Theory of Markoff Chains.  Transactions of the American Mathematical Society. 1942;52(1):37–64.
  * [64] Doob JL.  Markoff Chains–Denumerable Case.  Transactions of the American Mathematical Society. 1945;58(3):455–473.
  * [65] Gillespie D.  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  Journal of Computational Physics. 1976;22:403–434.
  * [66] Gillespie D.  Exact stochastic simulation of coupled chemical reactions.  The Journal of Physical Chemistry. 1977;81:2340–2361.  doi:10.1021/j100540a008.
  * [67] Redner S.  A guide to first-passage processes.  Cambridge University Press; 2001.
  * [68] Gardiner CW, et al.  Handbook of stochastic methods. vol. 3.  Springer Berlin; 1985.
  * [69] Balakrishnan V, Khantha M.  First passage time and escape time distributions for continuous time random walks.  Pramana. 1983;21(3):187–200.
  * [70] Khantha M, Balakrishnan V.  First passage time distributions for finite one-dimensional random walks.  Pramana. 1983;21(2):111–122.
  * [71] Phair RD, Scaffidi P, Elbi C, Vecerová J, Dey A, Ozato K, et al.  Global nature of dynamic protein-chromatin interactions in vivo: three-dimensional genome scanning and dynamic interaction networks of chromatin proteins.  Molecular and Cellular Biology. 2004;24(14):6393–6402.