{'fraction_non_alphanumeric': 0.0344023655180257, 'fraction_numerical': 0.00983737063573297, 'mean_word_length': 4.660444158352108, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phase Transition in a Model of Y-Shaped Molecules

In recent years the statistical mechanics of non-spherical molecules, such as polypeptide chains and protein molecules, has garnered considerable attention as their phase behavior has important scientific and health implications. One example is provided by immunoglobulin, which has a"Y"-shape. In this work, we determine the phase diagram of Y-shaped molecules on a hexagonal lattice through Monte Carlo Grand Canonical ensemble simulation, using histogram reweighting, multicanonical sampling, and finite-size scaling. We show that (as expected) this model is a member of the Ising universality class. For low temperatures, we implemented multicanonical sampling to induce faster phase transitions in the simulation. By studying several system sizes, we use finite-size scaling to determine the two phase coexistence curve, including the bulk critical temperature, critical chemical potential, and critical density.


I. INTRODUCTION
Our immune system is our primary defense against pathogenic organisms and malignant cells. A major component of this system is the class of proteins known collectively as immunoglobulin G (IgG). These are multidomain proteins that are particularly important due to their ability to bind to antigens with remarkable specificity. Antibodies have many practical applications due to this unusual specificity: this includes a large number of diagnostic applications in vitro, such as immunofluorescence, western blotting, and enzyme-linked immunosorbent assay analysis [1]. Antibodies have also become of intense interest in the pharmaceutical world, due to their efficacy as therapeutic molecules. This is evidenced by the large number of antibodies that are either approved, or in clinical trials for treating human disorders such as cancer, rheumatoid arthritis, osteoporosis, and asthma [2][3][4]. In such situations, highly concentrated antibody solutions are required in order to have therapeutic effect. However, such high concentration solutions present important stability and delivery challenges, such as aggregation and large solution viscosity. A large value of the viscosity occurs for IgG2, for example, at concentrations of 150 mg/ml, apparently due to a transition of the solution to a gel [5]. Therefore, understanding the condensation of proteins in solutions of highly concentrated immunoglobulin is of importance in pharmaceutical applications.
IgG molecules have a characteristic Y-shape, with two distal F ab arms that bind selectively to particular antigens. All of them have the same size and Y-conformation, due to their common genetic basis. However, they differ in their specific sequence of amino acids in the variable domains of the F ab arms. As noted by Wang et al. [6], the variability in the amino acid sequences can produce a large increase in the overall attractive interactions between neighboring IgG molecules. As a consequence, these interactions can cause a variety of phase transitions, including reversible aggregation, liquid-liquid phase separation, crystallization, and gelation, which occur in a wide variety of proteins [7][8][9]. Although IgGs are typically quite soluble at physiological conditions, sometimes they can become insoluble. In fact, recent studies of protein condensation have been published both for recombinant pharmaceutical IgGs and monoclonal IgGs [10][11][12][13][14][15][16][17]. A detailed discussion of the importance of IgGs in physiological and pharmaceutical situations is given by Wang et al. [6], as well as Nezlin [18]. Wang et al. stress that systematic studies of the phase behavior and phase diagrams of IgG solutions are crucial for the understanding of the pathological condensation in humans, as well as the stability of antibody drug formulations.
Of particular interest to this paper are the experimental studies of liquid-liquid phase transitions by Benedek et al [6,11,17]. They have reported several important results in a recent study of liquid-liquid phase separation in eight human myeloma IgGs and two recombinant pharmaceutical human IgGs. The first thing to note is that all liquid-liquid coexistence curves have quite similar shapes. These curves differ from those found for quasi-spherical proteins in that they are broad and asymmetric and have relatively small critical concentrations. The similarity of their shapes is presumably due to their common Y-shape. The critical temperatures of these liquid-liquid phase separations vary from one IgG to another, due to the variability in the amino acid sequences that leads to different net attractive protein-protein interactions. These transitions are also metastable and are preempted by the stable fluid-solid(freezing) transition, with the transition temperatures typically in the range of −20 • C to −30 • C. The study of such phase transitions is a relatively new, but rapidly emerging research area. It is therefore of interest to examine the effects of such an unusual architecture on model studies of phase transitions. One such study has already been carried out on a model of IgG [19].
In this paper we present an initial study of Y-shaped molecules by focusing on such molecules on two-dimensional lattices. Although this model might lack some important features of the phase transitions undergone by IgG in three dimensions (for example, it does not yield any notable asymmetry in the phase diagram), it does provide an example of a phase transition of a molecule with an unusual architecture. In addition, there are examples of IgG on two dimensional surfaces. For example, human IgG has been adsorbed to biomaterial surfaces, which can enhance long-term macrophage adhesion in vitro [20]. (However , in that case the Y molecules are not aligned only in the plane, but have some alignment perpendicular to the plane, in contrast to our model.) Due to the geometry of the molecule, we choose to use hexagonal lattices on which the Y-molecules can be naturally placed.
To study the possible phase transition of our model, we use grand-canonical Monte Carlo simulations along with the usual histogram reweighting and multicanonical biasing methods. The outline of our paper is as follows: In Section II we define our model and the method by which we study it. In Section III we present our results, namely that this model belongs to the Ising universality class. This particular result was expected a priori since the order parameter is a scalar. Using finite-size scaling (FSS) methods, we also obtain the liquid-liquid phase separation curve, including the location of the critical point in the thermodynamic limit. We also provide results for the finite-size system critical point, using a standard scaling analysis. In Section IV we present a brief conclusion.

II. MODEL AND SIMULATION
Our model of the Y-shaped molecule is implemented within the environment of a two-dimensional hexagonal lattice with side length L containing L 2 sites. Periodic boundary conditions are employed. Each of the N molecules occupies 4 lattice sites; the four include one central site and 3 rigid distal arms. This leads to a maximum number density of the system ρ of 0.25. The interactions included in this study are those between distal arms, and their respective nearest-neighbor distal arms, as shown in Fig.1, each with interaction strength of J. For simplicity, there are no center-to-center, centerto-distal arm, or any lattice binding energies. The total energy of the system U is then the summation of all distal armdistal arm interactions. We execute grand-canonical Monte Carlo (GCMC) simulations. To analyze the data and obtain the phase diagram, we use the Bruce-Wilding finite-size scaling (FSS) techniques [21,22], along with histogram reweighting and multicanonical sampling methods [23], to compile the phase diagram of this system. The fact that the order parameter for this model (defined below) is a scalar suggests that this Y-molecule model belongs to the Ising universality class, as we show in Section III.
Assuming that our model belongs to the Ising universality class, the critical point of our system can be determined by matching the probability density function (PDF) of the order- ing operator M of our system with the universal distribution of the two-dimensional Ising class. The order parameter M for the fluid is given by [21,22] where u = U/N is the energy density, and s and r are system specific parameters to be determined later. Similar to the order parameter M , the energy-like parameter ε is given by The Ising universality class has two relevant scaling fields, namely h, the ordering scaling field, and τ , the thermal scaling field. For fluids in this universality class, τ and h are defined as where ω = J/kT , µ is the reduced chemical potential in units of kT , and the subscript c denotes the critical point. The parameters r and s determine the degree of mixing in the relative scaling fields as well as M and ε.
During a simulation in a system size of side L, at fixed values of µ and ω, we record the molecule number density ρ and the energy density u, from which we determine the joint probability density function P (ρ, u). The joint PDF, P (M, ε), for the rescaled variables, M and ε, is related to the joint distribution of density and energy such that P (M, ε) = (1 − sr)P (ρ, u). We focus mostly on the order parameter M is a scaling parameter such that P M (x) has unit variance. Therefore, the PDF P L (x) of our model must also matchP M (x) at the fixed point. It is only for large L that the numerically obtained PDF tends to the fixed point distributionP M (x).
The fixed point functionP M (x) for the two-dimensional Ising model has been determined from previous work [21]. For one system size, by varying T , µ, and s, and matching the numerically obtained P L (x) withP M (x), we can determine the T c (L), µ c (L), and ρ c (L) of our model at said system size. By repeating this process for multiple system sizes, the bulk critical temperature(T c ), chemical potential(µ c ), and density(ρ c ) can be extrapolated. We obtain the parameter r from the slope of the µ-ω coexistence line at criticality [22] as seen in Fig.2.
To avoid performing numerous simulations, we use the standard method of histogram reweighting. To analyze the data and obtain the phase diagram, we use the Bruce-Wilding FSS techniques outlined here [21,22,24,25]. In order to obtain the critical parameters of the infinite system, we performed GCMC simulations for systems with side lengths 30, 40, 50, and 60 with periodic boundary conditions. The observables recorded during the simulation were u and ρ, from which, P (ρ, u), P (ρ), P (M, ε), and others were calculated. For each temperature, chemical potential, and system size, the simulation ran for 5,000 -6,000 Monte Carlo steps (MCS), and for 15,000 -25,000 Monte Carlo steps for simulations implementing and not implementing biasing techniques respectively, before recording the density and energy of the system. Each MCS comprises N attempts to change the system either by a molecule translation, rotation, insertion, or removal. The changes that the GCMC attempts to make to the system are defined as follows: Translation: Attempt to move the center of a randomly selected molecule to one of the 6 nearest neighboring sites of the center.
Rotation: Attempt to rotate a randomly selected molecule either clockwise, or counter-clockwise about the molecule center.
Insertion: Attempt to place a molecule, with a random orientation, on a randomly chosen lattice site.
Removal: Attempt to remove a randomly selected molecule.
This algorithm is ergodic. The density and energy of the system are recorded a total of 150,000 -250,000 times over the length of the simulation, with 1,000 -2,000 and 2,000 -4,000 MCS in between each recording for systems with and without the preweighting function, respectively.

III. RESULTS
The method we use to find the bulk critical parameters is as follows. First, we plot P L (x) for varying system sizes and find best fits to the universal fixed pointP M (x) by varying T c (L), µ c (L), and s. A sample best fit of P L (x) to this fixed point is shown in Fig.3 and shows FIG. 2: Plot of reduced chemical potential µ versus the reduced interaction strength ω. The parameter r is the slope of this curve at the critical temperature and is found to be r = −6.111.
our study, our model belongs to the Ising universality class. Next, we use the FSS predictions [24] that and to determine the bulk values of the critical temperature and critical chemical potential. In these equations, T c and µ c are the actual bulk critical temperature and chemical potential, and T c (L) and µ c (L) are the apparent bulk critical temperature and chemical potential determined from matching P L (x) to the fixed universal distributionP M (x). θ is a correction to scaling exponent. We use the value θ = 1.35, as calculated by Barma and Fisher[26], which coincides with the conjecture of Nienhuis [27] for the two-dimensional Ising system. The next step is to record T c (L) and µ c (L) for each system size, plot T c (L) and µ c (L) versus L −(θ+1)/ν , and then extrapolate to the infinite system size for both T c (L) and µ c (L). We then record the extrapolated points as T c and µ c . The resulting graphs of T c (L) and µ c (L) versus L −(θ+1)/ν are given in Fig.4. FSS also predicts a similar correction to the critical density of the model, namely where ρ c is the bulk critical density, ρ c (L) is the apparent critical density of a system of size L at T c and µ c , and d = 2 is the dimensionality of the system. Similar to the cases of T c and µ c , ρ c is extrapolated from the plot in Fig.5 and is ρ c = 0.127 ± 0.002. The coexistence curve was determined through a series of GCMC simulations at varying temperatures below the critical region, implementing the histogram reweighting and multicanonical sampling techniques discussed previously. Coexistence between two phases at a temperature T is confirmed when the areas underneath the two peaks in the density distribution P (ρ) are equal. The peak densities are recorded and plotted on the phase diagram. Examples of some density distributions for varying T at coexistence are shown in Fig.6.
In temperature regions apart from the critical region, finitesize effects are expected to be negligible since the correlation length is much smaller than that of the size of the system. With this expectation, the density peaks of P (ρ) found in our finite systems at different temperatures below the critical region will still mimic that of an infinitely large system. A phase diagram of the infinite system is constructed by determining the positions of the peaks in P (ρ) for sub-critical temperatures. As an additional check on this estimate of the equilibrium densities, we also calculate the average density of each the two phases present by using the PDFs to calculate these statistical averages. Our results are shown in Fig.7, where the FSS estimates of T c and ρ c are also plotted. We fit the points of our phase diagram to a power law of the form [24] This fit is also included in Fig.7. There is a small asymmetry evident on the high density side of the phase diagram, but not as pronounced as the one found experimentally or numerically in 3D for IgG. We show the critical parameters determined through the FSS method in Table I. We also show in Fig.8

IV. CONCLUSION
We end this paper with a few comments. First, the major feature of our work is that we have determined in detail using finite-size scaling methods the phase diagram of a molecule with unusual architecture. We have chosen a very simple model for the interaction between these molecules; it would be straightforward to include other interactions in future work. We have shown that the model belongs to the Ising universality class as one might expect, since its order parameter is a scalar. It seems clear that future research will increasingly deal with unusual molecular architectures. We note in passing the fact that the molecule has a Y-shape has not resulted in an asymmetry of the type found for IgG on the high density side. There are several reasons for this. First, our model is two-dimensional. Second, we have not taken into account in the model the large excluded volume effect that characterizes IgG, which Bendek et al. [6] have argued is responsible for its asymmetry. Indeed, the three-dimensional model of IgG studied by Li et al. [19] does have an excluded volume effect, and the phase diagram shows a pronounced asymmetry similar to that observed in experiment. It should also be noted that the reason our model has an asymmetry is due to the absence of a particle-hole symmetry, such as that present in the Ising model.
One direction for future research is to improve our model in order to describe the absorption of IgG on surfaces [19]. The inclusion of surfaces would require a somewhat more complex simulation algorithm, but would have the benefit of producing richer phase behavior. Another research direction would involve a more chemistry based, coarse-grained model of IgG, such as the one proposed by Voth's group [28], to study not only its phase transitions, but its viscosity. Such models are especially useful in probing the electrostatic interaction between antibodies and its impact on ordering. Finally, we note that the rheological information (such as the viscosity) that could be obtained with a coarse-grained model is of considerable interest in the pharmaceutical world.