1 Introduction
Genetic information is transmitted through genome and exhibits a hierarchical structure in the nucleus. The cell nucleus and the genome are organized into spatially separated subcompartments [erdel2018formation]. For eukaryotes, genetic information is stored in DNA which is also integrated into chromosomes of a cell nucleus [erdel2018formation, lee2017new, solovei2013lbr, cremer2006chromosome]. In chromosome territory, chromatin fibers are formed when DNA molecule is wrapped around the histones, and form transcriptionally inactivated and condensed region form as heterochromatin and transcriptionally activated region form euchromatin. The position, structure and configuration of heterochromatin and euchromatin regions are closely related to the gene expression. For instance, it has been observed that the volume of cell nucleus is a main determinant of the overall landscape of chromatin folding [gursoy2014spatial, gursoy2017spatial, gursoy2017computational]. Different kinds of nuclear architecture can attribute to different distribution and configuration of heterochromatin and euchromatin regions [strom2017phase, cremer2001chromosome], with additional proteinmediated specific interactions among genomic elements [perez2020chromatix, sun2021high]. Moreover, the inverted architecture occurs where heterochromatin is located at the center of the nucleus and euchromatin is enriched at the periphery. In contrast, when the heterochromatin is enriched at the nuclear periphery and around nucleoli, this is referred as the conventional architecture.
In [solovei2009nuclear, solovei2013lbr], Solovei et al. demonstrated that different types of nuclear architecture were associated with different mammalian lifestyles, such as diurnal versus nocturnal. The inverted nuclear architecture can be transformed from the convention architecture in mouse retina rod cells [solovei2009nuclear, solovei2013lbr]. The reorganization process is accompanied by the relocation of chromosomes from positions enriched at nuclear periphery, and the recreation of a single heterochromatin cluster into the inverted architecture. The difference of nuclear structure is partially attributed to the activity of lamin B, lamin A and envelope proteins [paulsen2017chrom3d, solovei2009nuclear, solovei2013lbr, kinney2018chromosome, van2017lamina]. Moreover, the rate of conversion of heterochromatin to euchromatin can also be controlled by volume constraints of the nucleus.
Recently several mathematical models using the phase field approach [laghmach2020mesoscale, lee2017new, seirin2017role] have been introduced to study the mechanism of generation of different nuclear architectures, including the size and shape of the nucleus, the rate of conversion of heterochromatin to euchromatin. These models usually include the minimization of total energy with various relevant geometric constraints. A common method to preserve such constraints is through extra penalty terms introduced in the energy functionals in models [laghmach2020mesoscale, lee2017new]. The drawback of such methods is the presence of the large penalty parameters which results in a stiff nuclear architecture reorganization systems, leading to significant challenge in the simulation and analysis. This is especially true in situations when the volume of chromosome must be enforced during the reduction of nuclear size and reorganization of nuclear architectures.
The Lagrange multiplier approach is commonly used for constrained gradient dynamic systems [du2006simulating, du2009energetic, du2008numerical, MR4049377, yang2017efficient, yang2021numerical, wang2016efficient, cheng2018multiple]. In this paper, we introduce a new Lagrange multiplier approach [cheng2020new, cheng2019global, sun2020structure] to enforce the geometric constraints such as the volume constraints for both chromosome and heterochromatin. When the volumes of each chromosome and heterochromatin are preserved as constants, the reactiondiffusion system with the Lagrange multipliers leads to a constrained gradient flow dynamics which satisfies an energy dissipation law. We develop several numerical schemes for nuclear architecture system with Lagrange multipliers. One is a weakly nonlinear scheme which preserves the volume constraints but requires solving a set of nonlinear algebraic systems for the Lagrange multipliers, the second is a purely linear scheme which approximates the volume constraints to second order and only requires solving linear systems with constant coefficient. These two schemes, while being numerically efficient, do not satisfy a discrete energy law. Hence, we construct the third scheme which is also weakly nonlinear but is unconditionally energy stable. However, this scheme requires solving a nonlinear algebraic system of (where is the number of chromosomes in the nucleus) equations which may require smaller time steps to be well posed. One can choose to use one of these schemes in different scenarios. In our numerical simulations, we use the first scheme to control the volumes to targeted values exactly, then we switch to the second scheme which is more efficient. We can use the third scheme if we want to make sure that the scheme is energy dissipative.
For validation purpose, we present several simulations results which are consistent with those observed in the experiments. We also demonstrate that our model and schemes are efficient and robust for investigating various nuclear architecture reorganization processes.
The paper is organized as follows: in Section 2, we present our phase field model for nuclear architecture reorganization by using a Lagrange multiplier approach. We chose suitable energy functionals to capture the most important interactions and constraints of various biological elements, and introduce Lagrange multipliers to capture the specific geometric constraints for the biological process. In Section 3, we develop efficient linear and weakly nonlinear time discretization schemes for the phase field model developed earlier. We present numerical results using the proposed schemes in Section 4, and compare them with the existing experimental literature and previous works. Finally we conclude the paper with more discussions of our methods and results in Section 5.
2 A phase field model for nuclear architecture reorganization (NAR)
To study the nuclear architecture reorganization (NAR) process, we employ a phase field/diffusive interface method. To start with this approach, the total nucleus is defined by a phase (labeling) function such that and respectively for the interior and exterior regions of the nucleus (Fig. 2). For a ellipsoid shape domain, we can define
(2.1) 
where is the interfacial (transitional domain) width, and describe the ellipse shape of the nucleus. Similarly, we will introduce phase functions to describe the heterochromatin region and to describe each individual chromosome region, where represents the total number of chromosomes in the nucleus ( Fig. 2). In particular, we will choose chromosomes for drosophila and chromosomes for human. In addition to the heterochromatin region, the rest of chromosome region is the euchromatin region. These are the order parameter/phase field functions that will be used to determined the final nuclear architecture.
For the general phase field approaches, the configuration and distribution of various regions are the consequence of minimizing a specific energy functional in terms of the above phase field functions, which takes into all considerations of the coupling and competition between different domains, as well as the relevant geometry constraints.
In the NAR models [lee2017new, laghmach2020mesoscale], the free energy for chromosome and heterochromatin had been chosen to include
(2.2) 
where and measure the interfacial thickness corresponding to and , is the computational domain, and is the double well potential which possess the local minima at and (see [provatas2011phase]). This part of free energy represents the competition and coupling between various chromosome and heterochromatin regions.
Next we will consider the following geometric constraints for all these regions that are biologically relevant to our application [alberts2015essential, cremer2010chromosome]:

All chromosomes are restricted within the cell nucleus region;

Heterochromosome of each chromosome stays within the chromosome;

Between all chromosomes, due to the excluded volume effects, do not selfcross or cross each other.
These constraints could be incorporated into the model by introducing three extra terms in the free energy:
(2.3) 
where and are three positive constants which indicate the intensities of domain territories, and is used for the induction of driving interface between and while keeping the local minima and fixed during dynamic process. The required conditions for are
(2.4) 
The lowest degree polynomials satisfying the above conditions is , see Fig. 1.
Due to the expression of LBR and lamin A/C in the nuclear envelope, interactions between heterochromatin and the nuclear envelope play an important role in nuclear architecture reorganization process. For this purpose, another term was introduced in the free energy to describe the interactions [solovei2013lbr]:
(2.5) 
where is the affinity constant, and implies the heterochromatin will locate at the nuclear periphery, while leads to the lack of heterochromatin and nuclear envelope interactions due to the absence of LBR or lamin A/C. More precisely, represents the intensity of affinity between nuclear function and heterochromatin region function .
2.1 NAR model with Lagrange multipliers
In the general nuclear reorganization process, often one need to take into account more geometric constraints. In particular, we shall include the following constraints in our model:

The nuclear space is fully occupied by chromosomes;

Heterochromatin is converted from/to euchromatin within each chromosome.
Notice that our approach could be extended to more general situations, especially related to those of item 4.
Our approach here is to introduce Lagrange multipliers to guarantee the constraints in the following NAR model with the total free energy
(2.6) 
The corresponding AllenCahn type gradient flow [liu2003phase, cheng2020new, feng2003numerical, chen1998applications, guan2014second, guan2017convergence, yang2017efficient] with respect to the above energy and the constraints 4 and 5 take the form:
(2.7)  
(2.8)  
(2.9)  
(2.10) 
where is mobility constant, and represent, respectively. From (2.9), the volume of chromosome can be contracted or expanded to a given volume by using our model. The volumes of th chromosome and heterochromatin in the th chromosome at time during nuclear reorganization (growth or inversion) stage are enforced by the Lagrange multipliers , (2.9). The boundary conditions can be one of the following two types
(2.11) 
where is the unit outward normal on the boundary .
Let and be, respectively, the target volumes for each chromosome and heterochromatin in each chromosome, can be interpreted as the heterochromatin conversion rate during nuclear architecture reorganization process. We assume that they will reach the target values at time and then stay there according to:
(2.12) 
where and are determined by and . and are suitable positive constants related to the time scale. Similarly, we assume that and evolve according to
(2.13) 
where and are determined by and with and being the targeted axis lengths for the ellipse enclosing the nucleus.
Remark 2.1
In [lee2017new], a penalty approach is introduced to satisfy the physical constraints 4 and 5 by adding the following to the free energy:
(2.14) 
where are three positive penalty parameters. The volume of th chromosome and volume of heterochromatin in th chromosome are defined in (2.9). A disadvantage of the penalty approach is large penalty parameters are needed for accurate approximation of the physical constraints, and may lead to very stiff systems that are difficult to solve numerically. The Lagrange multiplier approach that we use here can enforce these nonlocal constraints exactly and is free of penalty parameters. Furthermore, the NAR model (2.7)(2.10) based on the Lagrange multiplier approach can control exactly the growth rate of volume for different compartments during the nuclear reorganization (growth or inversion).
Remark 2.2
In the phase field approaches, there are many ways to represent the volume of each individual domain. For instance, the volume of each chromosome denoted by and their corresponding heterochromatin domain volume could be computed by the integrals and . However this representation may have disadvantages in the minimizing procedure, especially for the penalty methods used in [laghmach2020mesoscale, lee2017new, seirin2017role], due to its linearity with respect to the phase functions. One way to overcome this is to use the polynomial function (see Fig. 1) for the computation of the volumes for different chromosome regions. Since is an increasing function with respect to in the interval with and . One can adapt and for the volume of th chromosome and heterochromatin in th chromosome.
Let be the inner product in , and be the associated norm in . The constrained NAR model (2.7)(2.10) with (2.12) can be interpreted as a gradient system which implies phase separations will happen for .
Theorem 1
Proof
Taking the inner product of (2.7) with , we obtain
(2.16) 
Taking the inner product of (2.8) with , we obtain
(2.17) 
We derive from (2.9) that
(2.18) 
Similarly, we obtain
(2.19) 
Summing up (2.16) for , combing it with (2.17)(2.19) and using equality
(2.20) 
we obtain the desired energy dissipation law (2.15).
3 Numerical Schemes
In this section, we construct several efficient time discretization schemes based on the Lagrange multiplier approach [cheng2020new, cheng2019global] for the phase field NAR model (2.7)(2.10). For the sake of simplicity, for any function , we denote , and .
We split the total energy into a quadratic part and the remaining part as follows:
where is
(3.21) 
Once the volume of nucleus is given, volumes of each chromosome territory can be set up accordingly so that the constraint (2.10) can be satisfied automatically.
3.1 A weakly nonlinear volume preserving scheme
Note that in the first stage, we need to increase the volumes and to the targeted values and according to (2.12), respectively. Hence, we shall first construct below a volume preserving scheme which allows us to achieve this goal. More precisely, a second order scheme based on the Lagrange multiplier approach is as follows:
(3.22)  
(3.23)  
(3.24)  
(3.25) 
Below we show how to efficiently solve the coupled scheme (3.22)(3.25). Writing
(3.26) 
in (3.22), collecting all terms without , with or with , we find that, for , can be determined from the following decoupled linear systems:
(3.27) 
(3.28) 
(3.29) 
Then, writing
(3.30) 
in (3.23), we find that and can be determined from the following decoupled linear systems:
(3.31) 
(3.32) 
We observe that the above systems are all linear Poissontype equation with constant coefficients so they can be efficiently solved.
3.2 A linear scheme
In practice, the scheme (3.22)(3.25) should be used if we want to exactly preserve the volume dynamics of chromosome and heterochromatin. A disadvantage of the scheme (3.22)(3.25) is that we need to solve a nonlinear algebraic system which may require small time steps. To accelerate the simulation, we construct below a linear scheme for system (2.7)(2.10) which is more efficient but only approximately preserve the volume dynamics.
To this end, we reformulate (2.7)(2.10) into the following equivalent system:
(3.33)  
(3.34)  
(3.35)  
(3.36) 
Note that the last two relations are obtained by taking the time derivative of and in (2.9).
A secondorder linear scheme for the new system (3.33)(3.36) is as follows:
(3.37)  
(3.38)  
(3.39)  
(3.40) 
The above coupled scheme can be solved in essentially the same fashion as the scheme (3.22)(3.25). In fact, setting
(3.41) 
in (3.37), we find that for , are still determined from (3.27) (3.29). Then, writing
(3.42) 
in (3.38), we find that and are also determined from (3.31) (3.32). Once we have obtained and , we plug (3.41)(3.42) into (3.39)(3.40) to obtain a linear algebraic system for that can be solved explicitly. In summary, the scheme (3.37)(3.40) can be efficiently implemented as follows.
3.3 A weakly nonlinear energy stable scheme
Note that the schemes (3.22)(3.25) and (3.37)(3.40) are not guaranteed to be energy dissipative. Below we modify the scheme (3.22)(3.25) slightly to construct a weakly nonlinear but energy stable scheme with essentially the same computational cost for when volumes of each chromosome and heterochromatin become constants.
The idea is to introduce another Lagrange multiplier to enforce the energy dissipation. To this end, we introduce another Lagrange multiplier and expand the system (3.33)(3.36) for as
(3.43)  
(3.44)  
(3.45)  
(3.46)  
(3.47)  
Remark 3.3
Since volumes of each chromosome and heterochromatin are constants for , we have for . This zero term is critical for constructing energy stable schemes.
Then, a secondorder energy stable scheme based on system (3.43)(3.47) can be constructed as follows.