Chromatin state switching in a polymer model with mark-conformation coupling

TL;DR

This paper models chromatin as a polymer with coupled magnetic and conformational states, revealing a first-order phase transition that explains discrete chromatin state switching during cell differentiation.

Contribution

It introduces a polymer-Potts model with mark-conformation coupling, demonstrating a first-order phase transition in chromatin-like systems using mean-field analysis.

Findings

Phase transition is first-order with a large magnetization jump.

Small changes in epigenetic factors can cause macroscopic chromatin state switches.

Model explains discrete chromatin state changes observed in biological processes.

Abstract

We investigate the phase transition properties of the polymer-Potts model, a chain composed of monomers with magnetic degrees of freedom, with the motivation to study the conformation and mark switching dynamics of chromatin. By the mean-field approximation, we find that the phase transition between the swollen-disordered state and the compact-ordered state is discrete; it is first-order as in the long-range Potts model, but with a significantly larger jump in magnetization (i.e., mark coherence) upon the ordering transition. The results imply how small changes in epigenetic writer concentrations can lead to a macroscopic switching of the chromatin state, suggesting a simple mechanism of discrete switching observed, for instance, in cell differentiation.

Figures

Equations

$$f(\rho ,\{x_i\})$$

$$\frac{f(\varphi ,\{x_i\})}{k_{\mathrm{B}}T}$$

$${\tilde{J}}^{(\mathrm{c})}\simeq \tilde{v}+2\sqrt{\ln 2}$$

$${\varphi }^{\ast }\simeq \sqrt{(\tilde{J}-\tilde{v})/2}\left(\mathrm{for}\tilde{J}>{\tilde{J}}^{(\mathrm{c})}\right)$$

$$2{\left(\Delta \varphi \right)}^2=\sum _{1\le i,j\le q}\left({\tilde{J}}_{ij}^{(\mathrm{c})}+\Delta {\tilde{J}}_{ij}\right)\left[S_i\left(\tilde{\mathbf{h}}\right)+\Delta x_i\right]\left[S_j\left(\tilde{\mathbf{h}}\right)+\Delta x_j\right]-\tilde{v}.$$

$$\sum _{1\le i,j\le q}{J}_{ij}^{(\mathrm{c})}{S}_i\left(\tilde{\mathbf{h}}\right)S_j\left(\tilde{\mathbf{h}}\right)=v.$$

$$\overline{f}=\min \{\underset{\varphi ,m}{\min }f(\varphi ,m),\underset{\eta ,m}{\min }{f}^{\prime }(\eta ,m)\},$$

$$\frac{f^{\prime }(\eta ,m)}{k_{\mathrm{B}}T}$$

$$H=H_{\mathrm{GC}}(\{{\mathbf{R}}_n\})+H_{\mathrm{P}}(\{{\mathbf{R}}_n\},\{s_n\})+H_{\mathrm{M}}(\{s_n\})+H_{\mathrm{VE}}(\{{\mathbf{R}}_n\})+H_{\mathrm{F}}({\mathbf{R}}_0,{\mathbf{R}}_N).$$

$$H_{\mathrm{GC}}(\{{\mathbf{R}}_n\})=\frac{3k_{\mathrm{B}}T}{2b^2}\sum _{n=1}^N{\left({\mathbf{R}}_n-{\mathbf{R}}_{n-1}\right)}^2,$$

$$H_{\mathrm{P}}(\{{\mathbf{R}}_n\},\{s_n\})=-\sum _{n\ne m}J(s_n,s_m)\delta \left({\mathbf{R}}_n-{\mathbf{R}}_m\right),$$

$$H_{\mathrm{M}}(\{s_n\})=-\sum _i{h}_i\sum _n{\delta }_{s_n,i},$$

$$H_{\mathrm{VE}}(\{{\mathbf{R}}_n\})=v\sum _{n\ne m}\delta \left({\mathbf{R}}_n-{\mathbf{R}}_m\right)+w\sum _{n\ne m\ne k}\delta \left({\mathbf{R}}_n-{\mathbf{R}}_m\right)\delta \left({\mathbf{R}}_m-{\mathbf{R}}_k\right),$$

$$H_{\mathrm{F}}({\mathbf{R}}_0,{\mathbf{R}}_N):=-\mathbf{F}\cdot ({\mathbf{R}}_N-{\mathbf{R}}_0),$$

$$H_{\mathrm{P}}(\{{\mathbf{R}}_n\},\{s_n\})=-\sum _{i,j}{J}_{ij}\int \mathrm{d}\mathbf{r}{\sigma }_i(\mathbf{r},\{{\mathbf{R}}_n\},\{s_n\}){\sigma }_j(\mathbf{r},\{{\mathbf{R}}_n\},\{s_n\}),$$

$$H_{\mathrm{VE}}(\{{\mathbf{R}}_n\})=v\int \mathrm{d}\mathbf{r}\rho (\mathbf{r},\{{\mathbf{R}}_n\}{)}^2+w\int \mathrm{d}\mathbf{r}\rho (\mathbf{r},\{{\mathbf{R}}_n\}{)}^3,$$

$${\rho }_{\mathrm{av}}(\mathbf{r},{\mathbf{R}}_{N0}):=\frac{3N}{4\pi {R_{N0}}^3}\theta (R_{N0}-r),$$

$${\sigma }_{i,\mathrm{av}}(\mathbf{r},{\mathbf{R}}_{N0},\{s_n\}):={\rho }_{\mathrm{av}}(\mathbf{r},{\mathbf{R}}_{N0})\frac{1}{N}\sum _n{\delta }_{i,s_n},$$

$$H^{\mathrm{LR}}=H_{\mathrm{GC}}(\{{\mathbf{R}}_n\})+H_{\mathrm{P}}^{\mathrm{LR}}({\mathbf{R}}_{N0},\{s_n\})+H_{\mathrm{M}}(\{s_n\})+H_{\mathrm{VE}}^{\mathrm{LR}}({\mathbf{R}}_{N0})+H_{\mathrm{F}}({\mathbf{R}}_0,{\mathbf{R}}_N).$$

$$H_{\mathrm{P}}^{\mathrm{LR}}({\mathbf{R}}_{N0},\{s_n\})$$

$$H_{\mathrm{VE}}^{\mathrm{LR}}({\mathbf{R}}_{N0})$$

$$Z:=\left(\prod _{n=1}^N\int \mathrm{d}{\mathbf{R}}_n\sum _{s_n=1}^q\right)\exp \left(-\frac{H^{\mathrm{LR}}}{k_{\mathrm{B}}T}\right).$$

$$Z$$

$$\left(\prod _{n=1}^N\int \mathrm{d}{\mathbf{R}}_n\right)\delta (\mathbf{R}-{\mathbf{R}}_{N0})\exp \left[-\frac{3}{2b^2}\sum _{n=1}^N({\mathbf{R}}_n-{\mathbf{R}}_{n-1}{)}^2\right]$$

$$=\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\exp \left(\mathrm{i}\mathbf{k}\cdot \mathbf{R}\right)\left(\prod _{n=1}^N\int \mathrm{d}{\overline{\mathbf{R}}}_n\right)\exp \left(-\frac{3}{2b^2}\sum _{n=1}^N{{\overline{\mathbf{R}}}_n}^2-\mathrm{i}\mathbf{k}\cdot \sum _{n=1}^N{\overline{\mathbf{R}}}_n\right)$$

$$={\left(\frac{2\pi b^2}{3}\right)}^{3N/2}\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\exp \left(-\frac{N{b}^2{\mathbf{k}}^2}{6}+\mathrm{i}\mathbf{k}\cdot \mathbf{R}\right)$$

$$=\frac{1}{N^{3/2}}{\left(\frac{2\pi b^2}{3}\right)}^{3(N-1)/2}\exp \left(-\frac{3R^2}{2N{b}^2}\right).$$

$$Z$$

$$\left(\prod _{n=1}^N\sum _{s_n=1}^q\right)\exp \left(\frac{1}{R^3}\sum _{i,j}\frac{J_{ij}^{\prime }}{k_{\mathrm{B}}T}\sum _{n,m}{\delta }_{i,s_n}{\delta }_{j,s_m}+\sum _{i=1}^q\frac{h_i}{k_{\mathrm{B}}T}\sum _{n=1}^N{\delta }_{i,s_n}\right)$$

$$=\left(\prod _{i=1}^q\sum _{N_i=0}^N\right)\exp \left(\frac{1}{R^3}\sum _{i,j}\frac{J_{ij}^{\prime }}{k_{\mathrm{B}}T}{N}_i{N}_j+\sum _{i=1}^q\frac{h_i}{k_{\mathrm{B}}T}{N}_i\right)\left(\prod _{n=1}^N\sum _{s_n=1}^q\right)\left[\prod _{i=1}^q\delta \left(N_i,\sum _{n=1}^N{\delta }_{i,s_n}\right)\right]$$

$$=\left(\prod _{i=1}^q\sum _{N_i=0}^N\right)\exp \left(\frac{1}{R^3}\sum _{i,j}\frac{J_{ij}^{\prime }}{k_{\mathrm{B}}T}{N}_i{N}_j+\sum _{i=1}^q\frac{h_i}{k_{\mathrm{B}}T}{N}_i\right)\frac{N!}{{\prod }_{i=1}^q{N}_i!}\delta \left(\sum _{i=1}^q{N}_i,N\right).$$

$$\frac{N!}{{\prod }_{i=1}^q{N}_i!}=\exp \left(-\sum _{i=1}^q{N}_i\ln \frac{N_i}{N}-\frac{1}{2}(q-1)\ln N+O(N^0)\right).$$

$$\left(\prod _{i=1}^q\sum _{N_i=0}^N\right)\exp \left(\frac{1}{R^3}\sum _{i,j}\frac{J_{ij}^{\prime }}{k_{\mathrm{B}}T}{N}_i{N}_j+\sum _{i=1}^q\frac{h_i}{k_{\mathrm{B}}T}{N}_i\right)\frac{N!}{{\prod }_{i=1}^q{N}_i!}\delta \left(\sum _{i=1}^q{N}_i,N\right)$$

Full text

††preprint: APS/123-QED

Introduction.— Chromatin is a large polymer composed of monomers called nucleosomes, which are histone protein complexes wrapped with DNA Cortini et al. (2016); Fierz and Poirier (2019). The switching of cell states is encoded in the changes in epigenetics Waddington (1957) such as in the molecular and structrual changes in the chromatin. The chromatin states have been traditionally categorized into two: euchromatin and heterochromatin, which correspond to active (open) and inactive (closed) parts of the chromatin in terms of gene expression and accessibility. Consistent with this, recent chromosome conformation capture (Hi-C) experiments Lieberman-Aiden et al. (2009); Dixon et al. (2015) have identified the existence of two major compartments in the genome. The regions within the same compartment share similar marks (i.e., chemical modifications) in the nucleosomes and tend to interact with each other more frequently than across Di Pierro et al. (2017). Experiments have shown that different marks on the nucleosome induce distinct interactions due to the natural attraction and repulsion between nuclesomes Funke et al. (2016); Gibson et al. (2019), or by mediating proteins such as HP1 Bannister et al. (2001); Machida et al. (2018), and the polycomb repressive complexes Angel et al. (2011); Boettiger et al. (2016).

It has also been established that cell differentiation is accompanied by a large (megabase) scale transition in the compartments as well as changes in the states of epigenetic marks Dixon et al. (2015). The mechanism behind this switching, however, remains elusive. Previous modeling studies have assumed mark-dependent interactions between nucleosomes in order to explain the observed contact maps and 3D structures of the chromatin Jost et al. (2014); Barbieri et al. (2012). Other models have considered how the interaction between the marked histones lead to bistability in the coherent epigenetic marks Dodd et al. (2007); Sneppen et al. (2008); Jost (2014). A natural question is then how the interplay of chromatin chain dynamics and the kinetics of nucleosome modifications can lead to the drastic switching of compartments observed in differentiation.

To model chromatin polymer dynamics under stochastic modifications of nucleosomes, the polymer-Potts model has been considered Garel et al. (1999); Michieletto et al. (2016); Coli et al. . In this model, the random motion of the polymer chain is accompanied by monomer-monomer interactions that depend on the histone marks, and the histone marks can stochastically switch due to enzymatic reactions and histone turnover Deal et al. (2010); Dion et al. (2007). It has been numerically shown Michieletto et al. (2016, 2017a) that even for the Ising-type model, where there is essentially only two distinct states of the histone marks, there is a first-order-like transition between the swollen-disordered state, which corresponds to a loose polymer with spatially random marks, and the compact-ordered state, where the conformation is globular and the marks are coherent. This abrupt transition is likely due to the coupling between the conformation change and the epigenetic switchings, although a concrete theory is still lacking.

In this paper, we investigate the phase transition properties of the polymer-Potts model by considering a polymer chain in continuum space with stochastic histone mark exchange. Employing the Flory-type mean-field approximation for the dynamics of the chain, we write the pseudo free energy of the generic polymer-Potts model as a function of the order parameters representing the magnetization (i.e., mark coherence) and the polymer conformation. For the Ising-type interaction, the transition between the swollen-disordered state and a compact-ordered state is first-order, consistent with simulations and theories investigating mark dynamics on self-avoiding random walks Garel et al. (1999); Coli et al. . In the general case with multiple types of marks, we find that the jump in magnetization at the transition point is always larger in the polymer-Potts model compared with the Potts-model counterpart, and also obtain a criteria for the absense of a continuous transition. We further study the switching transition upon stretching of the chain, which serves as a simple model of force-induced epigenetic modification.

Model.— We consider a polymer model with Potts-like interactions between monomers. The interactions are mediated in a histone-mark-dependent way by proteins that we call readers, and the marks can change stochastically due to enzymatic reactions caused by the writers (Fig. 1). We assume that there are $q$ ($\geq 2$) types of marks. Extending the Flory-type mean-field approximation de Gennes (1975) to the present situation, the pseudo free energy (per monomer) at a temperature $T$ reads

[TABLE]

Here $\rho$ is the average monomer concentration given as $\rho\sim N/R^{3}$ with the polymer end-to-end distance $R$ and the total number of monomers $N$. The other variables, $\{x_{i}\}$, represent the $i$-th type mark occupation ($\sum_{i=1}^{q}x_{i}=1$). The first two terms in Eq. (1) correspond to the volume exclusion effect ($v,w>0$). Note that the parameters $v$ and $w$ represent the second and third virial coefficients, respectively. The third term originates from the Potts-like two-body interactions between monomers mediated by the readers. The detail of the interactions between the different types of marks is coded in $J_{ij}$, which is a real symmetric matrix. The fourth term is the entropy associated with the mark degrees of freedom. The last term represents the effect of external fields $\{h_{i}\}$, which describes how much a specific epigenetic mark is favored, reflecting, for example, the concentration of the histone modification enzymes, i.e., the writers. The equilibrium state is determined by minimizing Eq. (1) with respect to $\rho$ and $\{x_{i}\}$.

A few remarks are to be made for Eq. (1). Firstly, fluctuation effects are neglected compared with the microscopic model, although its inclusion will likely not change the key results Moore (1977); Lifshitz et al. (1978). Secondly, while higher-order interaction terms are irrelevant near the conventional second-order coil-globule transition point de Gennes (1975), inclusion of these terms may shift the transition point if the coil-globule transition becomes first-order, as in the situations explained below. Nevertheless, we expect that Eq. (1) captures the key characters observed in simulations of similar systems Michieletto et al. (2016); Coli et al. and is useful to generally analyse models with multiple kinds of marks. Lastly, the chromatin state transitions of our interest is at the level of sub-regions of compartments or several topologically associated domains Dixon et al. (2012), which is megabase scale corresponding to $N=10^{3-5}$. Although we have omitted all the terms that vanish in the limit of $N\to\infty$ in Eq. (1), it is straightforward to include higher-order terms and discuss their effects on the properties of transition de Gennes (1975); SM .

For the sake of understanding, let us first fix the mark degrees of freedom, $\{x_{i}\}$. Then, Eq. (1) is equivalent to the free energy of a classic polymer de Gennes (1975) in the large $N$ limit. As investigated in de Gennes (1975), there exists a transition between the coiled state, a swollen polymer with the average length scaling as $R\sim N^{3/5}$ Doi and Edwards (1986); de Gennes (1979), and the globule state, a densely packed polymer with $R\sim N^{1/3}$ de Gennes (1975), upon changing the overall two-body interaction (in the present case, $v-\sum_{i,j}J_{ij}x_{i}x_{j}$) from repulsive to attractive. Such coil-globule transitions have been observed in experiments using DNA Yoshikawa et al. (1996) and chromatin Zinchenko et al. (2017). The coil-globule transition in this case is continuous in the limit of $N\to\infty$ de Gennes (1975), even beyond the mean-field approximation Moore (1977); Lifshitz et al. (1978).

In another direction of simplification, we can consider the order-disorder transition of the marks under a fixed polymer density, $\rho$. Assuming a globular configuration ($R\sim N^{1/3}$) and $J_{ij}=J\delta_{ij}$, Eq. (1) represents the mean-field free energy of the Potts model Wu (1982). In the limit of $N\to\infty$, the order-disorder transition is continuous for $q=2$ while discontinuous for $q\geq 3$, which has been believed to be correct in three dimensions even beyond the mean-field approximation Wu (1982).

Phase transition by interactions and fields.— Introducing the dimensionless globular order parameter $\phi:=(w/k_{\mathrm{B}}T)^{1/4}\sqrt{\rho}$, we can express the equilibrium free energy as $\overline{f}=[\min_{\phi,\{x_{i}\}}f(\phi,\{x_{i}\})\text{ s.t.}\sum_{i=1}^{q}x_{i}=1]$, where

[TABLE]

Here, the dimensionless two-body interaction strengths are defined as $\tilde{J}_{ij}:={J}_{ij}/\sqrt{k_{\mathrm{B}}Tw}$ and $\tilde{v}:=v/\sqrt{k_{\mathrm{B}}Tw}$. As we have seen, if $\{x_{i}\}$ or $\phi$ is fixed to some value, the system described by Eq. (2) will show the conventional coil-globule or magnetic transition, respectively.

To see the effect of the coupling between $\{x_{i}\}$ and $\phi$, we first study the Ising-model case, $q=2$ and $J_{ij}=J(2\delta_{ij}-1)$ with $h_{i}=0$. Introducing the magnetization, $m:=x_{1}-x_{2}$, $f(\phi,m)$ for the case of $\tilde{v}=0.1$ is shown in Fig. 2. Let us denote the equilibrium values of $\phi$ and $m$ as $\phi^{*}$ and $m^{*}$, respectively. In this model, there is a critical value $\tilde{J}^{\mathrm{(c)}}$ such that for $\tilde{J}<\tilde{J}^{\mathrm{(c)}}$, the swollen-disordered phase is the equilibrium [$\phi^{*}=0$ and $m^{*}=0$, Fig. 2(a)], whereas for $\tilde{J}>\tilde{J}^{\mathrm{(c)}}$, this switches to the compact-ordered state [$\phi^{*}>0$ and $|m^{*}|>0$, Fig. 2(c)]. At the transition point [$\tilde{J}=\tilde{J}^{\mathrm{(c)}}$, see Fig. 2(b)], both the swollen-disordered state and the compact-ordered state are stable, meaning that there is a first-order transition. Thus, a switching transition can occur by simply changing the strength of the reader-mediated interaction, as has been seen numerically in a similar model Michieletto et al. (2016).

In Fig. 3(a) the $\tilde{J}$ dependence of $\phi^{*}$ and $m^{*}$ is plotted for $\tilde{v}=0.1$, showing a clear jump of the order parameters at the transition point. Since $|m^{*}|\simeq 1$ in the compact-ordered state (Fig. 2), we can approximate $\min_{\phi,m}f(\phi,m)\simeq\min_{\phi}f(\phi,m=\pm 1)=-(\tilde{J}-\tilde{v})^{2}k_{\mathrm{B}}T/4$. We then obtain

[TABLE]

shown as the black dotted line in Fig. 3(a), giving a good approximation. We have confirmed that the features such as the jump of $|m^{*}|$ from [math] to $\simeq\!\!1$ at the transition point are observed for a broad range of the values of $\tilde{v}$ SM .

We further consider the effect of external mark-specific fields by setting $h_{1}=h$ and $h_{2}=-h$ in the Ising-type model. The phase diagram in the case of $\tilde{v}=0.1$ is shown in Fig. 3(b). We find that the field-induced transition from the swollen-disordered state to the compact-ordered state is of first-order around the zero-field transition point, meaning that increasing or decreasing specific writers can also induce the switching behavior. Interestingly, within a certain range of $\tilde{J}$, sequential second- and first-order transitions occur as the field becomes stronger [Fig. 3(c)]. For a smaller $\tilde{J}$, a single continuous transition is induced by the field. Note that such a field-induced transition has been discussed in the context of the magnetic polymer within mean-field approaches Garel et al. (1999) as well as in simulations on self-avoiding walk models Garel et al. (1999); Huang and Luo (2004).

Propeties of transitions under general settings.— Here we consider the condition for the transition between the swollen and compact states to be countinuous under general $q$, $\{J_{ij}\}$, and $\{h_{i}\}$. Minimizing Eq. (2) on the assumption that the continuous transition occurs at $\{\tilde{J}_{ij}^{\mathrm{(c)}}\}$ for a given set of $\{h_{i}\}$, we obtain $\phi^{*}=0$ and $x_{i}^{*}=S_{i}(\tilde{\bm{h}}):=\exp(h_{i}/k_{\mathrm{B}}T)/\sum_{j=1}^{q}\exp(h_{j}/k_{\mathrm{B}}T)$ SM . The order parameters will grow in response to the deviation of the interaction strengths from their critical values: $\phi^{*}=\Delta\phi$ and $x_{i}^{*}=S_{i}(\tilde{\bm{h}})+\Delta x_{i}$ for $\tilde{J}_{ij}=\tilde{J}_{ij}^{\mathrm{(c)}}+\Delta\tilde{J}_{ij}$. Minimizing Eq. (2) at $\tilde{J}_{ij}=\tilde{J}_{ij}^{\mathrm{(c)}}+\Delta\tilde{J}_{ij}$ will give the following relation:

[TABLE]

For the transition to be continuous, the order parameters should smoothly change at the transition point: $\Delta\phi\to 0$ and $\Delta x_{i}\to 0$ for $\Delta\tilde{J}_{ij}\to 0$, meaning that $J_{ij}^{\mathrm{(c)}}$ should obey

[TABLE]

Therefore, if $\sum_{i,j}J_{ij}S_{i}(\tilde{\bm{h}})S_{j}(\tilde{\bm{h}})<v$ is satisfied, which is when the mean effective two-body interaction is repulsive, any continuous swollen-compact transition is prohibited and only switch-like transitions can occur. This condition can in principle be checked in experiment by measuring the effective interactions between nucleosomes.

A simple example is again the Ising-type model without external fields. In this case, since $\sum_{i,j}J_{ij}S_{i}(\tilde{\bm{h}})S_{j}(\tilde{\bm{h}})=\sum_{i,j}J_{ij}/4=0<v$, the continuous conformation transition is always prohibited, consistent with our numerical results that the first-order transition occurs irrespective of the value of $v$ SM .

To investigate the order-parameter jump at the first-order transition more specifically, we consider the Potts-type interaction: $J_{ij}=J\delta_{ij}$ with $h_{i}=0$. Figure 3(d) shows the magnetization jump $\Delta m=\max\{x_{i}^{*}\}-\min\{x_{i}^{*}\}$ for $q\geq 2$ in the mean-field polymer-Potts model [Eq. (2)] compared with the conventional mean-field Potts model. In the Potts model Wu (1982), the transition is of second-order for $q=2$ while first-order for $q\geq 3$, and $\Delta m$ is given as $(q-2)/(q-1)$. In the polymer-Potts model, on the other hand, the magnetic transition is always first-order, and $\Delta m$ is a monotonically increasing function of $v$ and $q$, while the dependency on $v$ is almost negligible for $q\geq 4$. Notice that $\Delta m$ in the polymer-Potts model is always larger than that in the corresponding Potts model, and $\Delta m$ is practically unity for $q\geq 4$. This suggests that the polymer conformation change that accompanies the magnetic transition reinforces the all-or-none switching property.

Mechanical discontinuous transition.— We here consider what happens when a stretching force is applied to the edges of a polymer chain with mark degrees of freedom. For simplicity, let us investigate the effects of an external force term added to the pseudo free energy [Eq. (2)] with the Ising-type interaction. The force term can be written as $f_{\mathrm{F}}=-\bm{F}\cdot\bm{R}/N$ with an external force $\bm{F}$ and the polymer end-to-end vector $\bm{R}$.

Within the mean-field level, the free energy including the effects of an external force is given as

[TABLE]

where $f^{\prime}(\eta,m)$ is another pseudo free energy including the effect of the external force:

[TABLE]

Here, the dimensionless polymer length, $\eta:=R/Nb$, can be interpreted as an order parameter characterizing a stretched state with the scaling $R\sim N$ Lai (1995). In Eq. (8), the second term represents the entropic elasticity de Gennes (1979), which is essential under stretched conditions. We denote the equilibrium point as ($\phi^{*},\eta^{*},m^{*}$). Since the globule state ($\phi^{*}>0$, $R\sim N^{1/3}$) and the stretched state ($\eta^{*}>0$, $R\sim N$) are incompatible, only one of ($\phi^{*},\eta^{*}$) can be finite and the other should be zero.

Figure 3(e) shows the changes in order parameters upon varying of the external force $F$ for the case of $\tilde{v}=0.1$ and $\tilde{J}=3$, in which the compact-ordered state ($\phi^{*}>0$, $\eta^{*}=0$, and $|m^{*}|>0$) is stabilized when $F=0$. We can see that a first-order transition occurs at a certain critical value $F^{\mathrm{(c)}}$, above which a stretched-disordered state ($\phi^{*}=0$, $\eta^{*}>0$, and $m^{*}=0$) emerges. The numerically obtained phase diagram in the $J$-$F$ plane is shown in Fig. 3(f). Note that the force-induced coil-globule transitions are believed to be discontinuous also in classical polymer models at $N\to\infty$ Halperin and Zhulina (1991); Lai (1995); Grassberger and Hsu (2002); Geissler and Shakhnovich (2002). In the polymer-Potts model, we find that the mark degrees of freedom become immediately disordered accompanying this stretching transition. Similar discontinuous transitions between a compact-ordered state and a stretched-disordered state have recently been seen in molecular dynamics simulations with short-range interactions Michieletto et al. (2017b).

Relation to molecular dynamics simulations.— To compare our results with the molecular dynamics simulations Michieletto et al. (2016) using Lennard-Jones-type interactions [Fig. 4(a)], we consider the virial expansion. By neglecting $\mathcal{O}(\phi^{6})$ terms and the existence of the neutral mark, we obtain the pseudo free energy [Eq. (2)] with additional terms proportional to $m^{2}\phi^{4}$ SM .

By minimizing the pseudo free energy, we obtain the optimized $\phi^{*}$ and $m^{*}$ as a function of the interaction strength $\epsilon/k_{\mathrm{B}}T$. Figure 4(b) shows that a discontinuous transition with a large jump of magnetization occurs at the interaction strength $\epsilon/k_{\mathrm{B}}T\simeq 0.85$, which is close to the simulation result Michieletto et al. (2016) ($\epsilon/k_{\mathrm{B}}T\simeq 0.9$ for $N=2000$). Although the virial expansion is not generally justfied for cases with a first-order transition, this result suggests that a simplified framework can connect the molecular level measurement of histone interactions Funke et al. (2016) to the compartment level chromatin state transition.

Discussion and conclusion.— Here we have studied the polymer-Potts model at the mean-field level and found that switch-like transitions are largely enhanced, compared with the transitions in a polymer model with unchangeable marks or the conventional Potts model, due to the coordination of the coil-globule and magnetic transition. The bistable property leading to the first-order transition fits with the phenomenology of chromatin state transition and cell differentiation. For instance, it has been shown that elimination of small kilobase-scale genome regions can induce compartment switching of a whole megabase-scale region Sima et al. (2019). This can be explained by the bistability of the chromatin state, which allows localized histone mark biases induced by transcription factors to spread macroscopically. The hysteresis effect, which is expected to accompany the chromatin discontinuous transition, may also improve the stability of the epigenetic regulation against chemical and mechanical perturbations and cell division.

Additional to the equilibrium phase transition scenario proposed in this paper, nonequilibrium features of the chemical dynamics Michieletto et al. (2016); Coli et al. and the phase separation properties of the key components in chromatin dynamics Gibson et al. (2019); Strom et al. (2017); Larson et al. (2017); Plys et al. (2018); Tatavosian et al. (2019) may play roles in enhancing or diminishing the switch-like behavior. Nevertheless, the fact that a simple mark-conformation coupling can lead to a discrete switch indicates that nonlinear dynamics and well-designed chemical networks may not be essential in explaining cell fate dynamics. In real differentiation, state switching occurs in sub-regions and does not expand to the whole chromosome Dixon et al. (2015). It is interesting to explore how specific regions in the genome set boundaries to prevent the phase transition dynamics from spreading into undesired regions Obersriebnig et al. (2016).

Acknowledgements.— We are grateful to Yohsuke T. Fukai and Soya Shinkai for fruitful discussions and the reading of the manuscript. This work was supported by JSPS KAKENHI Grant Numbers JP18H04760, JP18K13515, JP19H05275, JP19H05795.

I polymer-Potts model

In this section, we describe how to obtain the pseudo free energy of the polymer-Potts model starting a generic form of microscopic interactions. The Hamiltonian of the polymer-Potts model with external field and force terms reads

[TABLE]

Here, the degrees of freedom are the marks of each monomer, $s_{n}$ $(1\leq n\leq N)$, which represent the types of chemical modifications attached to the $n$-th histone along the chain, and the positions of the monomers, $\bm{R}_{n}$. We assume that there are $q$ $(\geq 2)$ types of histone marks that are distinguishable in the sense of interactions: $s_{n}\in\{1,2,...,q\}$. Also, we take $\{\bm{R}_{n}\}_{n=1}^{N}$ as a set of variables. The first term of Eq. (S1),

[TABLE]

corresponds to the Gaussian-chain interaction that restricts the positions of monomers adjacent to each other in the polymer, where $b$ corresponds to the Kuhn length (the size of the monomer), and $T$ is the temperature. The second term,

[TABLE]

describes the attractive or repulsive interaction between the monomers when they become close in contact with each other. The detail of interactions between the different types of marks can be coded in $J(s_{n},s_{m})=\sum_{1\leq i,j\leq q}J_{ij}\delta_{i,s_{n}}\delta_{j,s_{m}}$, where $J_{ij}$ is a real symmetric matrix. For example, the Ising-type interaction can be realized if we set $q=2$ and $J_{ij}=J(2\delta_{ij}-1)$, and the $q$-state Potts-type interaction can be realized if we choose $J_{ij}=J\delta_{ij}$. The third term,

[TABLE]

describes the effect of an external bias field, which favors a certain type of histone marks over other types. The fourth term,

[TABLE]

describes the two- and three-body volume exclusion interaction terms that we assume have no histone-mark dependence ($v,w>0$). The last term,

[TABLE]

represents a mechanical stretching force that is applied to the edges of the polymer.

I.1 Modified polymer-Potts model with long-range interactions

To perform a mean-field analysis of Eq. (S1), let us replace the local quantities with the averaged counterparts. First, we rewrite the interaction terms ($H_{\mathrm{P}}$ and $H_{\mathrm{VE}}$) with the polymer density, $\rho(\bm{r},\{\bm{R}_{n}\}):=\sum_{n}\delta(\bm{r}-\bm{R}_{n})$, and the mark density of the $i$-th state, $\sigma_{i}(\bm{r},\{\bm{R}_{n}\},\{s_{n}\}):=\sum_{n}\delta_{i,s_{n}}\delta(\bm{r}-\bm{R}_{n})$. We then obtain

[TABLE]

and

[TABLE]

by neglecting the terms proportional to $\delta(\bm{0})$. Next, we define the average density and mark density as

[TABLE]

and

[TABLE]

where $\bm{R}_{N0}:=\bm{R}_{N}-\bm{R}_{0}$ is the end-to-end distance vector and $\theta(x)$ is the Heaviside step function. Note that $\int\mathrm{d}\bm{r}\rho_{\mathrm{av}}(\bm{r},\bm{R}_{N0})=\int\mathrm{d}\bm{r}\rho(\bm{r},\{\bm{R}_{n}\})=N$ and $\int\mathrm{d}\bm{r}\sigma_{i,\mathrm{av}}(\bm{r},\bm{R}_{N0},\{s_{n}\})=\int\mathrm{d}\bm{r}\sigma_{i}(\bm{r},\{\bm{R}_{n}\},\{s_{n}\})=\sum_{n}\delta_{i,s_{n}}$.

As a mean-field approximation, we replace $\rho$ and $\sigma_{i}$ with $\rho_{\mathrm{av}}$ and $\sigma_{i,\mathrm{av}}$, respectively. The obtained Hamiltonian is a modified polymer-Potts model with long-range interactions:

[TABLE]

The long-range interaction terms ($H_{\mathrm{P}}^{\mathrm{LR}}$ and $H_{\mathrm{VE}}^{\mathrm{LR}}$) are given as

[TABLE]

where $J^{\prime}_{ij}:=3J_{ij}/4\pi$, and

[TABLE]

where $v^{\prime}:=3v/4\pi$ and $w^{\prime}:=(3/4\pi)^{2}w$.

I.2 Derivation of pseudo free energy

In the following, we derive the functional form of the pseudo free energy (per monomer) through analyzing the equilibrium free energy of the modified polymer-Potts model, $\overline{f}:=-(k_{\mathrm{B}}T/N)\ln Z$, where the canonical partition function is defined as

[TABLE]

Using the explicit form of $H^{\mathrm{LR}}$ and an identity, $\int\mathrm{d}\bm{R}\,\delta(\bm{R}-\bm{R}_{N0})=1$, we can obtain

[TABLE]

Note that $\bm{R}$ represents the end-to-end distance vector, which is clear from its definition.

To perform the integration with respect to $\{\bm{R}_{n}\}_{n=1}^{N}$, we change the integration variables from $\{\bm{R}_{n}\}_{n=1}^{N}$ to $\{\overline{\bm{R}}_{n}:=\bm{R}_{n}-\bm{R}_{n-1}\}_{n=1}^{N}$. Using $\delta(\bm{R}-\bm{R}_{N0})=\delta(\bm{R}-\sum_{n=1}^{N}\overline{\bm{R}}_{n})=(2\pi)^{-3}\int\mathrm{d}\bm{k}\exp[\mathrm{i}\bm{k}\cdot(\bm{R}-\sum_{n=1}^{N}\overline{\bm{R}}_{n})]$, we can proceed the calculation as

[TABLE]

We then obtain

[TABLE]

To perform the summation with respect to $\{s_{n}\}_{n=1}^{N}$, we use the identity, $\sum_{N_{i}=0}^{N}\delta(N_{i},\sum_{n=1}^{N}\delta_{i,s_{n}})=1$, where $\delta(a,b):=\delta_{a,b}$ is the Kronecker delta. We then obtain

[TABLE]

Note that $N_{i}$ represents the occupation number of the $i$-th mark state. Applying Stirling’s formula, we can estimate the terms originated from the entropy as

[TABLE]

Replacing the summation with respect to $N_{i}$ with an integration over the occupation ratio of the $i$-th mark state, $x_{i}:=N_{i}/N$, we obtain

[TABLE]

Based on the above calculation, after performing the solid-angle integration of $\bm{R}$, we finally obtain the following reduced form of the partition function:

[TABLE]

where the constant parts of the free energy, $f_{0}$, is defined as

[TABLE]

and the pseudo free energy, $f(R,\{x_{i}\})$, is given as

[TABLE]

In the limit of $N\to\infty$, the last term in Eq. (S23) is negligible, and we can apply Laplace’s method to estimate Eq. (S21). Thus, the free energy, $\overline{f}=-(k_{\mathrm{B}}T/N)\ln Z$, is obtained as

[TABLE]

Now we introduce two kinds of dimensionless order parameters, $\phi$ and $\eta$, regarding the conformation of the polymer. The first parameter, $\phi$, is defined as

[TABLE]

which corresponds to the square root of the globule density ($N/R^{3}$) and becomes finite in the compact (or globule) state satisfying $R\sim N^{1/3}$. The second parameter, $\eta$, is defined as

[TABLE]

which represents how the polymer is stretched and becomes finite in the stretched state satisfying $R\sim N$. Note that either of $\phi$ or $\eta$ can be finite in the limit of $N\to\infty$, which can be seen from the definitions of these order parameters. In the following, the order parameters that minimize the pseudo free energy (i.e., equilibrium state) are denoted as $\phi^{*}$ and $\eta^{*}$. Considering the scaling $R\sim N^{\nu}$, the swollen state ($\nu\simeq 3/5$ de Gennes (1979); Doi and Edwards (1986)) is included in the phase characterized by $\phi^{*}=\eta^{*}=0$. Although our framework cannot determine $\nu$ beyond the inequality $1/3<\nu<1$ within the swollen state in the $N\to\infty$ limit, $\nu=3/5$ is obtained in the parameter regime far from the transition point if we consider the finite size effect using Eq. (S23) de Gennes (1975).

Regarding the macroscopic mark, the occupation ratios of each mark, $\{x_{i}=N_{i}/N\}_{i=1}^{q}$, are order parameters. We denote the equilibrium mark order parameter as $x_{i}^{*}$, and call that the $i$-th mark is ordered when $x_{i}^{*}>1/q$.

Using $\phi$, $\eta$, and $\{x_{i}\}$, we can rewrite the formula of the free energy [Eq. (S24)] as

[TABLE]

The compact-state free energy, $\overline{f}_{\mathrm{C}}$, is defined as

[TABLE]

where the corresponding pseudo free energy is given as

[TABLE]

Here, $\tilde{v}:=v^{\prime}/\sqrt{k_{\mathrm{B}}Tw^{\prime}}=v/\sqrt{k_{\mathrm{B}}Tw}$ and $\tilde{J}_{ij}:=J^{\prime}_{ij}/\sqrt{k_{\mathrm{B}}Tw^{\prime}}=J_{ij}/\sqrt{k_{\mathrm{B}}Tw}$ are dimensionless interaction strengths, and

[TABLE]

is the non-interacting pseudo free energy of the marks with $\tilde{h}_{i}:=h_{i}/k_{\mathrm{B}}T$ being the dimensionless external field. On the other hand, the stretched-state free energy, $\overline{f}_{\mathrm{S}}$, is defined as

[TABLE]

where the corresponding pseudo free energy is given as

[TABLE]

Here, $\tilde{F}:=Fb/k_{\mathrm{B}}T$ is a dimensionless stretching force strength.

To sum up, in the case of $\overline{f}_{\mathrm{C}}<\overline{f}_{\mathrm{S}}$ with $\phi^{*}>0$ and $\eta^{*}=0$, the compact state ($R\sim N^{1/3}$) is stable and the realized order parameters minimize $f_{\mathrm{C}}(\phi,\{x_{i}\})$; in the case of $\overline{f}_{\mathrm{C}}>\overline{f}_{\mathrm{S}}$ with $\phi^{*}=0$ and $\eta^{*}>0$, the stretched state ($R\sim N$) is stable and the realized order parameters minimize $f_{\mathrm{S}}(\eta,\{x_{i}\})$; in the case of $\overline{f}_{\mathrm{C}}=\overline{f}_{\mathrm{S}}=\min_{\{x_{i}\}}f_{\mathrm{mark}}(\{x_{i}\})$ with $\phi^{*}=\eta^{*}=0$, the swollen state ($R\sim N^{3/5}$) is stable. The stable polymer state is also classified into ordered and disordered states according to the occupancy of each mark: the $i$-th mark is ordered (disordered) in the case of $x_{i}^{*}>1/q$ ($x_{i}^{*}\leq 1/q$).

I.3 Derivation of equations for the order parameters

Based on (S27), we derive the equations for the order parameters.

I.3.1 Compact state

To obtain the order parameters in the compact (or globule) state ($\phi^{*}>0$ and $\eta^{*}=0$), we need to minimize $f_{\mathrm{C}}(\phi,\{x_{i}\})$ subject to the constraint $\sum_{i=1}^{q}x_{i}=1$. Instead of directly treating the constraint, we introduce the Lagrange multiplier $\mu$. Then, as a necessary condition, the optimized order parameters should correspond to a stationary point of the following function:

[TABLE]

We can find the stationary points of $f^{\prime}_{\mathrm{C}}(\phi,\{x_{i}\},\mu)$ by solving $\partial f^{\prime}_{\mathrm{C}}/\partial\phi=0$, $\{\partial f^{\prime}_{\mathrm{C}}/\partial x_{i}=0\}_{i=1}^{q}$, and $\partial f^{\prime}_{\mathrm{C}}/\partial\mu=0$, i.e.,

[TABLE]

where $\phi^{*}$, $x_{i}^{*}$, and $\eta^{*}$ are the solution of these equations.

Since we here consider the compact state, we look for a solution with $\phi^{*}>0$. Then, Eq. (S34) reduces to

[TABLE]

where we use matrix representations, $[\bm{x}]_{i}:=x_{i}$ and $[\tilde{\mathsf{J}}]_{ij}:=\tilde{J}_{ij}$. On the other hand, Eq. (S35) can be rewritten as $x_{i}^{*}=\mathrm{e}^{\,\mu^{*}-1}\exp\,(2{\phi^{*}}^{2}\sum_{j=1}^{q}\tilde{J}_{ij}x_{j}^{*}+\tilde{h}_{i})$. Substituting this expression for $x_{i}^{*}$ in Eq. (S36), we obtain $\mathrm{e}^{\,\mu^{*}-1}=[\sum_{i=1}^{q}\exp\,(2{\phi^{*}}^{2}\sum_{j=1}^{q}\tilde{J}_{ij}x_{j}^{*}+\tilde{h}_{i})]^{-1}$. Therefore, we obtain the following expression of $x_{i}^{*}$:

[TABLE]

where $[\tilde{\bm{h}}]_{i}:=\tilde{h}_{i}$, and the softmax function, $S_{i}\,(\bm{y})$, is defined as

[TABLE]

As a consequence, the order parameters in the compact state, $\phi^{*}$ and $\{x_{i}^{*}\}_{i=1}^{q}$, satisfy Eqs. (S37) and (S38).

I.3.2 Stretched state

In the stretched state ($\phi^{*}=0$ and $\eta^{*}>0$), we need to minimize $f_{\mathrm{S}}(\eta,\{x_{i}\})$ subject to $\sum_{i=1}^{q}x_{i}=1$. Noticing that $\eta$ and $\{x_{i}\}$ are uncoupled, we can easily obtain

[TABLE]

and

[TABLE]

where $S_{i}\,(\bm{z})$ is the softmax function defined in Eq. (S39).

I.3.3 Swollen state

In the swollen state ($\phi^{*}=\eta^{*}=0$), we only have to minimize $f_{\mathrm{mark}}(\{x_{i}\})$ subject to $\sum_{i=1}^{q}x_{i}=1$. The solution of the order parameters, $\{x_{i}^{*}\}_{i=1}^{q}$, is already obtained in Eq. (S41).

We can derive a sufficient condition for the swollen state to be stabilized as explained in the following. Assuming $F=0$, we can see $\eta^{*}=0$ since from Eq. (S32), $f_{\mathrm{S}}(\eta,\{x_{i}\})$ is a monotonously increasing function of $\eta$. We thus show below a sufficient condition for $\phi^{*}=0$ on the assumption that $F=0$. Focusing on the stationary condition of $\phi^{*}$ [Eq. (S34)], we can see that $\phi^{*}$ has to be zero if $\sum_{i,j}\tilde{J}_{ij}x_{i}x_{j}={\bm{x}}^{\mathrm{T}}\,\tilde{\mathsf{J}}\,\bm{x}$ is smaller than $\tilde{v}$ for arbitrary values of $\{x_{i}\}_{i=1}^{q}$ satisfying $0\leq x_{i}\leq 1$. Introducing an orthogonal matrix $M$ which diagonalizes the symmetric matrix $\tilde{\mathsf{J}}$ as $\mathsf{M}^{\mathrm{T}}\,\tilde{\mathsf{J}}\,\mathsf{M}=\mathsf{\Lambda}:=\mathrm{diag}(\lambda_{1},\cdots,\lambda_{q})$, writing the maximum eigenvalue of $\mathsf{\Lambda}$ as $\lambda_{\mathrm{max}}$, and defining $\bm{y}:=\mathsf{M}^{\mathrm{T}}\bm{x}$, we obtain the following inequality,

[TABLE]

The second inequality is due to $0\leq x_{i}\leq 1$. Therefore, if the maximum eigenvalue of $\tilde{\mathsf{J}}$, $\lambda_{\mathrm{max}}$, is smaller than $\tilde{v}$, the only solution of Eq. (S34) is $\phi^{*}=0$. In short, a sufficient condition for the swollen state to be realized is that the maximum eigenvalue of $\mathsf{J}$ is smaller than $v$ (i.e., interaction is effectively repulsive) and also $F=0$ (i.e., no stretching force).

I.4 Criteria for second-order transition

Let us consider the case with some external field ($\bm{h}\neq\bm{0}$) and no external force ($F=0$), and explore under what conditions the second-order transition between a swollen state ($\phi^{*}=0$) and a compact state ($\phi^{*}>0$) can occur.

Recall first that in the swollen state, $\phi^{*}=0$ and $x_{i}^{*}=S_{i}(\tilde{\bm{h}})$ (Sec. I.3.3). Then, if we assume that a continuous transition occurs at $\tilde{J}_{ij}=\tilde{J}_{ij}^{\mathrm{(c)}}$, the system will smoothly change from the swollen state to the compact state in response to a certain small change in $\tilde{J}_{ij}$: $\phi^{*}=0+\Delta\phi$ and $x_{i}^{*}=S_{i}(\tilde{\bm{h}})+\Delta x_{i}$ for $\tilde{J}_{ij}=\tilde{J}_{ij}^{\mathrm{(c)}}+\Delta\tilde{J}_{ij}$. The pseudo free energy for $\tilde{J}_{ij}=\tilde{J}_{ij}^{\mathrm{(c)}}+\Delta\tilde{J}_{ij}$ is obtained from Eq.  (S29) as

[TABLE]

Minimizing this function with respect to $\Delta\phi$, we obtain

[TABLE]

The assumption of a continuous transition leads to the continuity of the order parameters at $\tilde{J}_{ij}=\tilde{J}_{ij}^{\mathrm{(c)}}$: $\Delta\phi\to 0$ and $\Delta x_{i}\to 0$ for $\Delta\tilde{J}_{ij}\to 0$; therefore, we obtain the following expression of the possible second-order transition point:

[TABLE]

This means that within the modified polymer-Potts model [Eq. (S11)], any second-order conformation transition accompanying the changes in the order parameters should satisfy Eq. (S45). Conversely, if Eq. (S45) is not satisfied, the continuous swollen-compact phase transition cannot occur. We will see in the following subsection I.5 that a discontinuous transition can occur even if Eq. (S45) is not satisfied.

I.5 Analysis of simple cases

In this subsection, we provide some examples to explain how the equilibrium states and transition properties of the polymer-Potts model can change according to the histone-mark specificity of the interactions, $\{J_{ij}\}$. For simplicity, we assume $F=0$ and $h_{i}=0$.

I.5.1 Mark-independent interaction

Let us consider interactions which are not specific to histone marks ($J_{ij}=J$). In this case, the formula of the free energy [Eq. (S27)] is reduced to $\lim_{N\to\infty}\overline{f}=\overline{f}_{\mathrm{C}}+\lim_{N\to\infty}f_{0}$ with the following pseudo free energy:

[TABLE]

The third term originated from the entropy of the marks can be minimized with $x_{i}=1/q$, so that we obtain the further reduced expression, $\lim_{N\to\infty}\overline{f}=\min_{\phi}f_{\mathrm{C,U}}(\phi)+\mathrm{const.}$, where the pseudo free energy for a uniform interaction, $f_{\mathrm{C,U}}(\phi)$, is defined as

[TABLE]

Minimizing this function, we can see that there is a continuous transition as the interaction strength $\tilde{J}$ is varied through a critical point, $\tilde{J}^{\mathrm{(c)}}=\tilde{v}$, and obtain the stable equilibrium states as follows:

[TABLE]

The obtained value of $\tilde{J}^{\mathrm{(c)}}$ is consistent with the discussion in the subsection I.4 [take $\bm{h}=\bm{0}$ and $\tilde{J}_{ij}^{\mathrm{(c)}}=\tilde{J}^{\mathrm{(c)}}$ in Eq. (S45)]. Figure S1 shows the $\tilde{J}$ dependence of the order parameter $\phi^{*}$ for the case of $\tilde{v}=0.1$.

Note that the mean-field model of the coil-globule transition introduced in de Gennes (1975) is reduced to the same functional form as $f_{\mathrm{C,U}}(\phi)$ in the limit of $N\to\infty$ (see also Sec. II). Thus, with a uniform Potts-like interaction, the marks are always disordered ($x_{i}^{*}=1/q$), and the properties of the polymer-Potts model are reduced to those of a usual polymer model with two-body and three-body monomer-monomer interactions.

I.5.2 Mark-specific interaction: Ising-type

Here and in the main text, we consider one of the simplest histone mark-specific interaction, the Ising-type interaction: $J_{ij}=J(2\delta_{ij}-1)$ with $q=2$. According to the general discussion in the subsection I.4 [see Eq. (S45)], there is no second-order transition in this system since the possible second-order transition point cannot be achieved irrespective of the value of $\tilde{J}=J/\sqrt{k_{\mathrm{B}}Tw}$.

It is natural to redefine the mark order parameter with the magnetization, $m:=x_{1}-x_{2}$. The numerically obtained order parameters are shown in Fig. S2. In stark contrast to the non-mark-specific interaction case discussed in the subsection I.5.1, there is a first-order transition between the swollen-disordered state ($\phi^{*}=0$, $m^{*}=0$) and the compact-ordered state ($\phi^{*}>0$ and $|m^{*}|>0$) at a $\tilde{v}$-dependent transition point, $\tilde{J}^{\mathrm{(c)}}$.

Since the magnetization amplitude $|m^{*}|$ shows a jump from [math] to $\simeq\!\!1$, we can obtain an approximate formula of the first-order transition point $J^{\mathrm{(c)}}$ and that of the $\tilde{J}$ dependence of $\phi^{*}$ in the compact state, as explained in the following. According to Eqs. (S27), (S28), and (S29), the free energy of the swollen-disordered state is given as $\overline{f}_{\mathrm{C}}=f_{\mathrm{C}}(\phi=0,x_{1}=x_{2}=1/2)=-\ln 2$. On the other hand, the free energy of the compact-ordered state is given as $\overline{f}_{\mathrm{C}}\simeq\min_{\phi}f_{\mathrm{C}}(\phi,x_{1}=1,x_{2}=0)\,[=\min_{\phi}f_{\mathrm{C}}(\phi,x_{1}=0,x_{2}=1)]=\min_{\phi}[(\tilde{v}-\tilde{J})\phi^{2}+\phi^{4}]$. The last expression can be minimized with $\phi=\sqrt{(\tilde{J}-\tilde{v})/2}$, and finally we obtain in the compact-ordered state $\overline{f}_{\mathrm{C}}\simeq-(\tilde{J}-\tilde{v})^{2}/4$. At the transition point ($\tilde{J}=\tilde{J}^{\mathrm{(c)}}$), the free energy of the swollen-disordered state and that of the compact-ordered state coincide with each other: $-\ln 2=-(\tilde{J}^{\mathrm{(c)}}-\tilde{v})^{2}/4$, and we thus obtain $\tilde{J}^{\mathrm{(c)}}\simeq\tilde{v}+2\sqrt{\ln 2}$.

To summarize, the stable equilibrium state for the Ising-type interaction is given as follows:

[TABLE]

Let us shortly consider the case with finite fields ($h_{1}=h$ and $h_{2}=-h$). Under sufficiently high fields, the possible second-order transition point is given from Eq. (S45) as

[TABLE]

As $h$ gets larger and the mark type becomes virtually fixed, $\tilde{J}^{\mathrm{(c)}}(h)$ gets close to $\tilde{v}$ Garel et al. (1999), which is the second-order transition point of a homogeneous polymer with only mark-independent interactions [see Sec. I.5.1 and Eq. (S48)]. In addition, for low fields where the magnetization jump is expected around unity as in the zero-field case, we can estimate the first-order transition point as in Eq. (S51):

[TABLE]

These expressions [Eqs. (S54) and (S55)] are plotted in Fig. 3(b) in the main text and compared with the numerical results.

Lastly, we consider external-force effects. In this case, the first-order transition occurs as a function of $F$ and $J$. Since the magnetization jump is almost unity, we can estimate the first-order transition point in a similar way to Eq. (S51):

[TABLE]

which is also compared with the numerical results [Fig. 3(f) in the main text].

I.5.3 Mark-specific interaction: $q$-state Potts-type interaction

As another example of a mark-specific interaction, let us think of the $q$-state Potts-type interaction: $J_{ij}=J\delta_{ij}$ with $q$ $(\geq 2)$ kinds of marks. In contrast to the Ising-type interaction case, according to Eq. (S45), there is a possible second-order transition at

[TABLE]

In Fig. S5, the $\tilde{J}$ dependence of the order parameters, $\phi^{*}$ and $m^{*}:=x^{*}_{1}-x^{*}_{2}$, is shown for the case of $q=2$. For small values of $\tilde{v}$ [Figs. S5(a) and (b)], two-step transitions occur: the compact-disordered state ($\phi^{*}>0$, $m^{*}=0$) starts to appear at $\tilde{J}=\tilde{J}^{\mathrm{(c)}}_{\mathrm{2nd}}$ through a second-order transition, and then the compact-ordered state ($\phi^{*}>0$ and $|m^{*}|>0$) emerges at $\tilde{J}=\tilde{J}^{\mathrm{(c)}}$ through a first-order transition. For large values of $\tilde{v}$ [Fig. S5(c)], on the other hand, a single-step transition occurs: the compact-ordered state emerges at $\tilde{J}=\tilde{J}^{\mathrm{(c)}}$ through a first-order transition. Defining the magnetization as $\Delta x^{*}:=\max_{i}\{x_{i}^{*}\}-\min_{i}\{x_{i}^{*}\}$, we can see that such qualitative features of the order parameters seem to be invariant even when the number of different marks ($q$) is larger than $2$ (see Figs. S5 and S5 for the case of $q=3$ and $q=10$, respectively).

For the case of large $\tilde{v}$, noticing that the magnetization size ($\Delta x^{*}$) shows a jump from [math] to $\simeq\!\!1$ we can estimate the first-order transition point ($\tilde{J}^{\mathrm{(c)}}$), just as in the case of the Ising-type interaction (see the subsection I.5.2). Comparing the free energy in the swollen-disordered state and that of the compact-ordered state, we can obtain

[TABLE]

On the other hand, for the case of small $\tilde{v}$, a similar comparison of the free-energy difference between the compact-disordered state and the compact-ordered state leads to $-(\tilde{J}/q-\tilde{v})^{2}/4-\ln q=-(\tilde{J}-\tilde{v})^{2}/4$, and we then obtain

[TABLE]

The approximate cross point ($\tilde{v}^{\mathrm{(cross)}}$, $\tilde{J}^{\mathrm{(cross)}}$) of the first-order and second-order transition lines can be estimated with Eqs. (S57) and (S58) [or Eqs. (S57) and (S59)] as

[TABLE]

The numerically acquired phase diagram in the $J$-$v$ plane for the case of $q=2$ is shown in Fig. S6 with the analytic expressions of the transition lines [Eqs. (S57), (S58), and (S59)].

II Polymer model with long-range interactions

In this section, we consider a polymer model with long-range interactions and show that this model is essentially equivalent to the mean-field model of the coil-globule transition introduced in de Gennes (1975).

Let us consider a polymer model with long-range interactions:

[TABLE]

where the Gaussian-chain term, $H_{\mathrm{GC}}(\{{\bm{R}}_{n}\})$, is defined in Eq. (S2), and the long-range volume exclusion interaction term, $H_{\mathrm{VE}}^{\mathrm{LR}}(\bm{R}_{N0})$, is defined in Eq. (S13). In this section, we assume that $v$ can be positive or negative, while $w$ is positive.

The equilibrium partition function, $Z_{\mathrm{polymer}}$, is defined as

[TABLE]

Introducing the end-to-end vector $\bm{R}$ with an identity, $\int\mathrm{d}\bm{R}\,\delta(\bm{R}-\bm{R}_{N0})=1$, we can rewrite $Z_{\mathrm{polymer}}$ in a similar way to Eq. (S15) as

[TABLE]

where $v^{\prime}=3v/4\pi$ and $w^{\prime}=(3/4\pi)^{2}w$ as in Sec. I. Through the same calculation process as Eq. (S16), we obtain

[TABLE]

Performing the solid-angle integration of $\bm{R}$ and exponentiating all the factors lead to the following form:

[TABLE]

Following de Gennes (1975), we define the expansion factor, $\alpha:=R/\sqrt{N}b$, which leads to

[TABLE]

where $F_{\mathrm{eff}}(\alpha)$ is an effective free energy defined as

[TABLE]

with $x:=v^{\prime}\sqrt{N}/k_{\mathrm{B}}Tb^{3}$ and $y:=2w^{\prime}/k_{\mathrm{B}}Tb^{6}$. The effective free energy, $F_{\mathrm{eff}}(\alpha)$, has the same functional form as the mean-field free energy, Eq. (1) in de Gennes (1975), apart from the numerical coefficient of the second term [$-3\ln\alpha$ in de Gennes (1975) instead of $-2\ln\alpha$ in Eq. (S68)]. This numerical difference of the coefficient only produces quantitative difference in the coil-globule transition properties for finite $N$, and moreover, the term proportional to $\ln\alpha$ becomes irrelevant in the limit of $N\to\infty$ since $\ln\alpha=o(N)$ (see also the discussion in the subsection I.5.1). Therefore, the polymer model with long-range two-body and three-body interactions is essentially equivalent to the mean-field free energy introduced in de Gennes (1975).

III Relation to microscopic models

In this section, we employ the virial expansion method to obtain a relation between the pseudo free energy of the polymer-Potts model and a microscopic Hamiltonian describing a polymer with mark-dependent interactions.

III.1 Gas model

Let us consider the microscopic two-body interaction term in the Hamiltonian of a virtual gas of $q$ types of monomers:

[TABLE]

Here, $\bm{r}_{n}$ and $s_{n}$ are the position and mark of the $n$-th monomer, respectively, and $U(\bm{r}-\bm{r}^{\prime};s,s^{\prime})$ is the interaction energy between one monomer with the mark $s$ placed at $\bm{r}$ and another monomer with the mark $s^{\prime}$ placed at $\bm{r}^{\prime}$. A basic inversion symmetry, $U(\bm{r}_{n}-\bm{r}_{m};s_{n},s_{m})=U(\bm{r}_{m}-\bm{r}_{n};s_{n},s_{m})=U(\bm{r}_{m}-\bm{r}_{n};s_{m},s_{n})$ is assumed. In the following, $U(\bm{r}_{n}-\bm{r}_{m};s_{n},s_{m})$ is also expressed as $U_{nm}$ for simplicity.

III.2 Virial expansion

Extending the standard procedure of the virial expansion Mayer and Mayer (1977) to the multiple-mark case, we can formally expand the interaction part of the partition function $Q:=V^{-N}\int\prod_{n=1}^{N}\mathrm{d}^{3}\bm{r}_{n}\exp(-H_{\mathrm{int}}/k_{\mathrm{B}}T)$ as

[TABLE]

where

[TABLE]

is the Mayer f-function.

Using an abbreviation, $U^{(ij)}_{12}:=U(\bm{r}_{1}-\bm{r}_{2};i,j)$ and correspondingly $f^{(ij)}_{12}:=f(\bm{r}_{1}-\bm{r}_{2};i,j)$, we can obtain the second-order virial coefficient as

[TABLE]

Using similar expressions, we can also obtain the third-order virial coefficient as

[TABLE]

Higher-order virial coefficients can be obtained in a similar manner.

III.3 Relation to polymer-Potts model with Ising-type interaction

Here, let us consider the Ising-type case where $q=2$ and $U(\bm{r};1,1)=U(\bm{r};2,2)$ (and thus $f_{12}^{(11)}=f_{12}^{(22)}$). Taking $v$ and $J$ as

[TABLE]

and

[TABLE]

we can readily obtain from Eq. (S72)

[TABLE]

where the density $\rho=N/V$ and the magnetization $m=(N_{1}-N_{2})/N$ are introduced. The form of Eq. (S76) corresponds to the two-body interaction terms in the polymer-Potts pseudo free energy [Eq. (1) in the main text] if we take $V\sim R^{3}$ with the polymer length $R$.

In a similar way, taking $w$ and $w^{\prime}$ as

[TABLE]

and

[TABLE]

we can obtain from Eq. (S73)

[TABLE]

The first term of Eq. (S79) corresponds to the three-body interaction term in the polymer-Potts pseudo free energy [Eq. (1) in the main text].

To sum up, the interaction part of the partition function can be represented as

[TABLE]

Thus, replacing $V=\mathrm{const.}\times R^{3}$ and adding the entropic contribution, we obtain the mean-field pseudo free energy (per monomer) of the polymer-Potts model in $N\to\infty$ limit:

[TABLE]

Here, the dimensionless globular order parameter $\phi=(w/k_{\mathrm{B}}T)^{1/4}\sqrt{\rho}$ and the dimensionless coupling constants $\tilde{v}=v/\sqrt{k_{\mathrm{B}}Tw}$, $\tilde{J}=J/\sqrt{k_{\mathrm{B}}Tw}$, and $\tilde{w}^{\prime}=w^{\prime}/w$ are used as in the main text. Note that the term proportional to $m^{2}\phi^{4}$ does not appear in the minimal model considered in the main text. Neglecting the higher-order terms $O(\phi^{6})$ (although, strictly speacking, this operation is not justified when a discontinuous coil-globule transition occurs), we can discuss based on Eq. (S81) the phase transition nature of the polymer-Potts model with general interactions satisfying the Ising symmetry [$U(\bm{r};1,1)=U(\bm{r};2,2)$].

III.4 Virial expansion of specific model

In the simulation of Michieletto et al. (2016), three kinds of marks are assumed ($1\leq i\leq 3$): two of them ($i=1$, $2$) with Ising-type interactions explained below and one of them ($i=3$) with no mark-specific interactions. For simplicity, we here neglect this third mark although its existence may shift the transition point slightly, as can the $O(\phi^{6})$ terms. The polymer model with two kinds of marks, where the interaction between a monomer with the mark $i$ and another monomer with the mark $j$ ($1\leq i,j\leq 2$) is given as the truncated Lennard-Jones type Michieletto et al. (2016):

[TABLE]

and $U_{\mathrm{LJ}}^{(ij)}(\bm{r})=0$ (for $r>r_{\mathrm{c}}^{(ij)}$). Here, $\sigma$ corresponds to the interaction range, $\epsilon_{ij}=\epsilon\delta_{ij}+k_{\mathrm{B}}T(1-\delta_{ij})$ represents the mark-specific interaction strength, $r_{\mathrm{c}}^{(ij)}=1.8\sigma\delta_{ij}+2^{1/6}\sigma(1-\delta_{ij})$ is the mark-specific interaction cutoff, and $\mathcal{N}=1+4(1.8^{-12}-1.8^{-6})$ is chosen so that $\min_{\bm{r}}U_{\mathrm{LJ}}^{(ii)}(\bm{r})=-\epsilon$ [see Fig. S7(a)].

Noticing the Ising symmetry of $U_{\mathrm{LJ}}^{(ij)}(\bm{r})$, i.e., $U_{\mathrm{LJ}}^{(11)}(\bm{r})=U_{\mathrm{LJ}}^{(22)}(\bm{r})$, we identify $U_{\mathrm{LJ}}^{(ij)}(\bm{r})$ as $U(\bm{r};i,j)$ in Sec. III.3 and calculate the dimensionless coupling constants $\tilde{v}$, $\tilde{J}$, and $\tilde{w}^{\prime}$ based on Eqs. (S74), (S75), (S77), and (S78). Since the Ising-type interaction strength can be quantified by a dimensionless parameter $\epsilon/k_{\mathrm{B}}T$, we plot $\tilde{v}$, $\tilde{J}$, and $\tilde{w}^{\prime}$ as a function of $\epsilon/k_{\mathrm{B}}T$ in Fig. S7(b). The mark-specific interaction strength $\tilde{J}$ is an increasing function of $\epsilon/k_{\mathrm{B}}T$ as expected, while $\tilde{v}$ and $\tilde{w}^{\prime}$ also depend on the value of $\epsilon/k_{\mathrm{B}}T$ since $\epsilon$ also affects the mark-nonspecific interactions.

By minimizing Eq. (S81), we can obtain the equilibrium globular order parameter $\phi^{*}$ and magnetization $m^{*}$, as shown in Fig. S7(c) as a function of $\epsilon/k_{\mathrm{B}}T$. There is a first-order transition between the swollen-disordered state ($\phi^{*}=0$, $m^{*}=0$) and the compact-ordered state ($\phi^{*}>0$, $|m^{*}|>0$) at $\epsilon/k_{\mathrm{B}}T\simeq 0.85$, which is of the same order as the critical value ($\epsilon/k_{\mathrm{B}}T\simeq 0.9$) obtained in the simulation with $N=2000$ Michieletto et al. (2016).

Lastly, we stress that the necessary condition for a second-order coil-globule transition derived in the main text [see the discussions around Eq. (6)] is applicable even when $m^{2}\phi^{4}$ terms or $O(\phi^{6})$ terms are present in the pseudo free energy as Eq. (S81). Indeed, $\tilde{v}$ is positive around the transition point [see Figs. S7(b) and (c)], and thus the second-order coil-globule transition is prohibited, consistent with the numerical result.