A modified Ising model of Barabási–Albert network with gene-type spins

The central question of systems biology is to understand how individual components of a biological system such as genes or proteins cooperate in emerging phenotypes resulting in the evolution of diseases. As living cells are open systems in quasi-steady state type equilibrium in continuous exchange with their environment, computational techniques that have been successfully applied in statistical thermodynamics to describe phase transitions may provide new insights to the emerging behavior of biological systems. Here we systematically evaluate the translation of computational techniques from solid-state physics to network models that closely resemble biological networks and develop specific translational rules to tackle problems unique to living systems. We focus on logic models exhibiting only two states in each network node. Motivated by the apparent asymmetry between biological states where an entity exhibits boolean states i.e. is active or inactive, we present an adaptation of symmetric Ising model towards an asymmetric one fitting to living systems here referred to as the modified Ising model with gene-type spins. We analyze phase transitions by Monte Carlo simulations and propose a mean-field solution of a modified Ising model of a network type that closely resembles a real-world network, the Barabási–Albert model of scale-free networks. We show that asymmetric Ising models show similarities to symmetric Ising models with the external field and undergoes a discontinuous phase transition of the first-order and exhibits hysteresis. The simulation setup presented herein can be directly used for any biological network connectivity dataset and is also applicable for other networks that exhibit similar states of activity. The method proposed here is a general statistical method to deal with non-linear large scale models arising in the context of biological systems and is scalable to any network size.


Introduction
Biological networks are multi-dimensional complex systems whose collective interaction in response to perturbations may lead to critical transitions from one stable state to another. Acute asthma attacks, clinical depression, diabetes mellitus, inflammation, and epileptic seizures, among many others, are examples of such sudden shifts in the state of the system, from healthy to diseased states (Trefois et al. 2015;Wolf et al. 2018). Such "phase transitions" are common in other complex systems such as ecological systems (for example, spontaneous extinction of species in response to gradual changes in external conditions) (Capitan and Cuesta 2010) or the evolution of human language (for example, the formation of Zipfian properties) (DeGiuli 2019; Vera et al. 2020). However, the idea that a system consisting of simple interacting units may exhibit phase transition was initially motivated by a seminal model of magnetic systems called the Ising model. The Ising model is one of the simplest and most frequently studied models of cooperative phenomena in statistical mechanics (Ising 1925). The classical Ising model is a discrete, pairwise interacting two-state system proposed to explain the structure and properties of ferromagnetic materials and has been solved exactly for one-and two-dimensional lattices (Onsager 1994). In the Ising model of a two-dimensional lattice, each site carries a spin which may be up or down, and neighboring spins prefer to be parallel to each other. The external field prefers to orient the spins in the direction of the field. The spins align in the same direction at low temperature, and the system exhibits spontaneous magnetization. At high temperatures, the spins align randomly, and the system is paramagnetic.
Since then, Ising models have been extended to study phase transitions occurring in more complicated topologies such as random, small-world and scale-free networks (Albert and Barabasi 2002;Bianconi 2002; Barrat and Weigt 2000;Dorogovtsev et al. 2002;Ferreira et al. 2010;Gitterman 2000;Herrero 2002Herrero , 2008Lopes et al. 2004;Pekalski 2001). For example, Ising models of networks can explain how the opinion of the individual is influenced by their contacts on opinions of people on a given subject (Aleksiejuk et al. 2002;Bartolozzi et al. 2006;Castellano et al. 2009;Contucci 2007;Redner 2017). Further real-world applications of Ising models of networks include socioeconomic problems such as racial segregation in the US, group herding, human culture (Stauffer and Hohnisch 2006); phase transitions in neural networks (Aldana and Larralde 2004); communication in the World Wide Web (Kumar et al. 2000); and systems biology (May and Lloyd 2001;Pastor-Satorras et al. 2015;Pastor-Satorras and Vespignani 2001).
In this regard, the analogy between phase transitions occurring in living systems (such as normal to diseased state transition) and physical systems (such as condensation of water) has been well-motivated (Davies et al. 2011;Holstein et al. 2013;Trefois et al. 2015;Smith 2010). The normal state to cancer state transition has been described as a process similar to the first-order irreversible discontinuous phase transition occurring in physical systems (Facciotti 2013;Jin et al. 2017;Liu et al. 2013;Mojtahedi et al. 2016;Torquato 2010). The central idea is that living systems are open systems in quasisteady state type equilibrium in continuous exchange with their environment wherein cells behave like a network in heat bath under external perturbations (Pastor-Satorras et al. 2015;Scheffer et al. 2012). They survive by exporting entropy to the environment in exchange for structural order. When a control parameter increases entropy, it causes collective flipping of states, which drives the system to an unstable critical state (or diseased state), thereby leading to phase transitions in living systems. In an Ising model, an analog for such a control parameter could be temperature or magnetic field, which, after a particular critical value, may cause the system to undergo a phase transition (discussed in detail in Sect. 2).
An integral part of such living systems is the biological networks that they are composed of-for example, the gene regulatory networks, protein interaction networks, among many others. A gene regulatory network represents a network of genes that can activate or suppress each other's expression levels owing to the interaction between the genes or due to the influence of agents external to the cell. One of the widely accepted method and a powerful tool for qualitative analysis of dynamics in gene networks is the Boolean dynamic modeling method. In a Boolean model of a network, nodes may be gene or protein and may either take on or off, indicating their expression levels, concentration, or activity. The relationships between the states of nodes are updated by logical functions or truth tables such as AND, OR, or NOT. Though this is a powerful tool, it requires that all possible states of the system i.e. 2 N (where N is the network size) be explicitly calculated for the time evolution of the network.
Further, most complex step scales as O(N 2 ) in the update schemes, and to our knowledge, the model is usually applied to networks whose sizes are of the order of hundred nodes (Campbell and Albert 2014;Wang et al. 2012;Zhang et al. 2008). This is because building such state transition graphs gets computationally expensive as network size increases due to the exponential dependence of the size of state space on network size, thereby making it challenging to analyze large-scale interconnected biological networks such as, for example, the complete human genome. As a consequence of this, when Boolean models are used for constructing signaling pathways on large and dense networks, the number of optimal solutions explodes, which necessitates alternative techniques from statistical physics and graph theory (Alexppoulos et al. 2010;Mitsos et al. 2009).
For large scale simulations, a simpler model that considers gene network as a system in continuous fluctuation that takes into account the current state of the nodes; and does not depend very much on the microscopic details or specific genomic features and is scalable to large sizes would be appropriate. In the method we propose, we overcome the issue with the scale by changing the way we update the state of nodes in response to external perturbations (such as temperature or magnetic field). Firstly, we consider the gene regulatory network as a system that exists in a quasi-steady state in thermal equilibrium with the environment, which aims to conserve energy as a whole whenever it changes its configuration. We establish an initial configuration for all nodes in the network, which may be all zeros or all ones or a combination of zeros and ones based on prior information. Then we perform a random walk over the configuration space. Specifically, this means that we randomly pick a node and calculate the energy cost for the system if this node were to change its state. The walker then hops from the initial configuration to this new configuration only if the transition probability is energetically favorable for the whole system, else the system retains its initial configuration ( (Metropolis et al. 1953) summarized in "Appendix"). Such "clever moves" are repeated multiple times to obtain averages to get the behavior of the system to the external perturbation. The behavior of the system is characterized by the mean of the summed states of the system, referred to as magnetization in the context of the classical Ising model. This method is an abstraction of the gene network and, therefore, only requires the initial configuration of all nodes of the network and network connectivity. It can calculate large scale response features, interpret gene expression of cells in large repositories, and is scalable to network size. This is a qualitative method that can give insights into overall organizing principles in the network and is capable of predicting coherently working clusters in the network. However, it is important to note that, in essence, the model we propose and the Boolean network model are not directly comparable since the former is a thermodynamic model, and the latter is a logical model. They are merely different ways to look at the same system. In the following Sect. 2, we introduce the modified Ising model and motivate the biological analogs of the Ising model by drawing parallels between physical and biological systems.

Model description and biological motivation
To describe phase transition in a system, we need to take into account interactions between its parts (Goldenfeld 1992;Kardar 2007). Since systems exhibit universal behavior near the critical point (Torabi and Davidsen 2019;Torabi and Rezaei 2016), a variety of statistical mechanical systems can be simulated by Ising-like models provided that the symmetry properties of the system, the pattern of interaction, and the dimensionality of the system is considered (Landau and Lifshitz 1980).
Given that the Ising model is simple and can predict cooperative behavior wherein each element has two states where the energy of each element depends on its state and that of its neighbors, it has found itself wide applications in addressing complexity in biology in the last century. The central aspect of these studies is that usually, the control parameters and physical properties of the Ising model are amenable to a biological interpretation depending on the target biological system being modeled.
In protein science, for example, a relatively popular adaptation of the Ising model, specifically the one-dimensional variant, is the homozipper, which is used as a simplified statistical thermodynamic model for protein folding. The homozipper is a sequence of multiple repeat proteins where each element of a repeat protein is an identical and independently folding unit that interacts with each other in a nearest-neighbor pairwise manner. The sequence can then be pulled apart like a zipper by mechanical force modeled by temperature. The folding is then a process constrained by the number of identical repeats, the energy of the repeated unit, and the interaction energy between the folded units given by the Hamiltonian of the Ising model (Aksel and Barrick 2009;Millership et al. 2016). This proposition of the Ising model to study order-disorder transitions in protein science extended from one-dimensions to higher dimensions for studying helix to coil transitions, beta-hairpin formation, hydrophobicity in protein chains and downhill folding (Garcia-Mira et al. 2002;Kubelka et al. 2004;Kubelka and Kubelka 2014;Lai et al. 2015;Naganathan and Munoz 2014;Munoz et al. 1997;Irback et al. 1996;Irback and Sandelin 2000;Lobanov and Galzitskaya 2011;Zimm and Bragg 1959).
In immunology, an Ising spin-model equivalent of the idiotypic-anti-idiotypic immunological networks has been shown to exhibit self-organization i.e. formation of large homogeneous domains at high temperatures. In such a system, each spin interacts with its mirror-image spin and the neighbors of the image. In the Ising model of such a system, spin up is synonymous to a proliferation of lymphocytes in an ocean of virgin states while spin-down represents a challenge to the immune system. Thus at low temperatures where there are few challenges and low noise, the system exhibits order. While at high-temperature i.e. lots of diseases, a disordered system is formed wherein the net magnetization synonymous to the activity level of lymphocytes is close to zero (Sahimi and Stauffer 1993).
The applications of the Ising model to study genome organization is not new. Ignoring the unique attributes of individual genes, (Baran and Ko 2006) has shown that transcription polarity in a bacterial chromosome i.e. the preference of genes to be coded in the leading strand of replication and their nature to form co-organized clusters can be modeled by a one-dimensional Ising model. Like the magnetic forces that align adjacent spins when the external magnetic field is applied, one could imagine adhesive pseudo-forces such as nearest-neighbor interaction that cause transcription orientation. The chromosome is then simply a series of spin-like indicator variables like that of the 2 N configurations that of the one-dimensional Ising lattice. Each gene is oriented negative or positive, depending on the sign of the open reading frame from which it is transcribed. Such models also allow themselves to be analytically tractable for the study of the effects of gene insertion and deletion. A one-dimensional longrange Ising model has been shown to be a rather robust description of long-range correlations in DNA sequences (Colliva et al. 2014).
Ising models have been used to analyze genetic data from affected sib-pair (ASP) where a data point can be represented as ±1 corresponding to an allele being shared or not by a sib-pair. The nearest neighbor interaction between adjacent dipoles is analogous to the interaction between adjacent genetic markers on a chromosome. The effect of an applied magnetic field i.e. a point field acting on a given particle is analogous to the effect of a disease gene, causing an increase in allele sharing at nearby locations. For example, in the ASP analysis, the coupling constant and magnetic field are interpreted as the strength of genetic linkage between markers and the effect of disease locus in distorting allele sharing in response to random genetic and environmental effects (Majewski et al. 2001).
Furthermore, Ising-and Potts-based models have been proposed for studying phase transitions occurring in more complicated topologies e.g., conformational restructuring affected by temperature. Here a single DNA strand is modeled as a system if interacting bases with short-and long-range interactions. Further, a set of such DNA strands live on a Cayley tree, which is a tree-like graph where each node has an equal number of branches. The edges of this graph may then take multiple spin values say, ±1, 0 where the former shows the existence of a Holliday junction, and zero means vacant or no edge. The response of the system to temperature is then analyzed based on a Translation Invariant Gibbs Measure (TIGM) (Rozikov 2017(Rozikov , 2018. To our knowledge, these studies do not find any direct mapping between the Ising parameters such as the temperature, magnetic field, or Boltzmann constant to genetic parameters. For an excellent review of similarities between physical and biological systems, we point the interested reader to (Davies et al. 2011).
More examples of applications of Ising models to biological systems include, but is not limited to, a four-dimensional cellular automaton-like Ising model in which cells transition between normal, proliferative, hypoxic and necrotic states has been used to model the tumorigenesis process which involves a transition between these pre-malignant and malignant cell states (Durrett 2013;Torquato 2010); estimating information transfer between spins occurring in human connectome (Marinazzo et al. 2014); the transition of B-DNA to S-DNA (Ahsan et al. 1998); estimation of differentially expressed genes in cancer patients (Xumeng et al. 2011); and approximation of join expression profiles of genes using a small number of observations (Santhanam et al. 2009).
However, to our knowledge, these models do not take into account two aspects of modeling phase transitions in biological networks that we address in this study. Firstly, the discovery of almost scale-free topologies of biological networks in the last couple of decades (Albert and Barabasi 2002). Secondly, the asymmetric i.e. 0, 1 states of activity that has provided a reasonable approximation of the reality of states in single cells (Cesar-Razquin et al. 2018;Wang et al. 2012). In this study, we address these aspects by proposing a modified Ising model for scale-free networks with genetype spins. The modified Ising model is an arrangement of genes or proteins in a cell where the interaction between the units could be short-range or long-range given by the connectivity matrix of the network. Each unit can interact with their connections in a pairwise manner. These units may activate or suppress each other, given by the binary state variable. A cell can survive in a healthy state if the majority of the units are in an active state, and if it gains energy from the external environment. However, external perturbations may slowly switch the cell to diseased states where the majority of the units may be inactive. Cell survival is, therefore, dependent on conserving energy and changing configurations only if it has a low cost. The Ising analogs for these external variables-energy and entropy are magnetic field and temperature, respectively. Broadly, these control parameters drive the phase transitions process occurring in the living systems.
In this paper, we establish a numerical and theoretical framework on a simulated scale-free network whose nodes exhibit binary states of activity. We show the conditions under which this network of gene-type spins undergoes phase transition due to the influence of temperature and magnetic field. This framework serves as a benchmark for future studies that aim to test dynamics of the Ising model of biological networks with gene-type spins from public databases. We refer to the model that comes out of this as the modified Ising model; a comparison of this model to the classical Ising model (Ising 1925) is summarized in Table 1.
Concretely, the Hamiltonian of the Ising model of such a network reads, where J is the coupling constant specifying the strength of interactions; A i j is the adjacency matrix; h indicates the constant external field; A i j s i s j is the coupling energy arising due to the interaction between nodes and shows the effect of cooperative behavior; h i s i is the energy arising due to the effect of magnetic field. The Hamiltonian so formed from these two terms is the total energy of the system. If J > 0, neighboring spins prefer to take the same values (referred to as ferromagnetic exchange interaction in a classical Ising model); when J < 0, neighboring spins prefer to take opposite values (referred to as anti-ferromagnetic exchange interaction in a classical Ising model). Spins, s i , s j can take values ±1 in the classical Ising model; and 0 and The higher h it is, the higher the interaction of the cell with the environment. This parameter (which corresponds to the magnetic field in the ferromagnetic Ising model) tries to retain active genes based on its interactions outside the cell. Minimizing energy in Eq. 3 corresponds to having active genes or healthy state as a result of those mentioned above pair-wise and environmental interactions. The temperature, on the other hand, induces fluctuation that results in randomness in the state of genes. The order parameter is then defined as the number of active genes in the cell, In the next sections, we will see how this system experiences a first-order phase transition due to the change in the parameter h whose critical value influence on the properties of the connectivity structure of the genes in the cell. Though the modified Ising model is proposed here with an eye on gene and protein interaction networks, the observations made here, in principle, hold for other similar real-world networks as well. Preliminary results of this work have been presented in the form of a poster and talk (Krishnan et al. 2018(Krishnan et al. , 2019Krishnan 2019). The paper is organized as follows: Sect. 2 provides a short overview of the Ising model and terminologies used in the subsequent sections of the paper along with biological motivation; in Sect. 3 we show the conditions under which the modified Ising model can undergo phase transitions for different initial configurations of the system (for positive and negative coupling constants) using Monte Carlo simulations; Sect. 4 presents the mean-field solutions and shows a mapping between classical Ising model of scale-free networks and the modified Ising model.

Numerical simulations
As motivated in Sect. 1 the focus of this paper is to study the system in Eq. 3 for modified Ising spins of the Barabási-Albert network. Such a network is constructed based on two main properties of a real-world network -linear growth and preferential attachment (Albert and Barabasi 2002). The network is initialized with m 0 nodes that are not connected. Subsequently, new nodes with m edges are added by an iterative growth process to the existing m 0 nodes. The resultant network has a power-law degree distribution and is characterized by a degree exponent, 2 < γ < 3 that resembles realworld biological networks.
We now put modified Ising spins on nodes of a Barabási-Albert network of size, N = 5 × 10 3 and preferentially attached links, m = 5. We choose this network size since it lies in a similar order of magnitude (≈ 10 3 ), such as that of gene regulatory network of standardized datasets such as S. cerevisiae or E. coli (Balaji et al. 2006;Gama-Castro et al. 2008). The number of links attached to grow the Barabási-Albert network is a free parameter and cannot, to our knowledge, be directly compared to real-world networks. Nevertheless we show the effect of the Barabási-Albert model parameters on the modified Ising model in the subsequent sections (cf. Fig. 8). Then with the standard heat bath Monte Carlo algorithm, we do a spin search for thermal equilibrium at temperature T (cf. algorithm in "Appendix"). The number of equilibration and sampling steps is proportional to the size of the network and has been chosen such that the system has had sufficient time to evolve from its initial configuration and reach a steady state. In other words, the system has visited different states in the phase space and can now generate states that are consistent with the parameters controlling the system. This can be verified by plotting properties of the system, such as magnetization, until it plateaus at a fixed value. We equilibrate the system for 2 × 10 4 MC steps and after this temporary period, we simulate for 3 × 10 4 MC steps. This allows for an average 10 spin flips per spin. We then sample at the end of every step and perform simulations for both ferromagnetically and antiferromagnetically coupled networks, under the influence and absence of the magnetic field.
Under no influence of the magnetic field and ferromagnetic exchange interaction, all nodes in the network start at an active state where the order parameter, M = 1. At T < 1, the system favors order as seen in the top panel of Fig. 1. As the thermal fluctuations in the system increases, the disorder in the system increases. The order parameter reaches 1 2 asymptotically as T → ∞. Similarly, when the system is initialized with an anti-ferromagnetic exchange interaction, the order parameter asymptotically reaches 1 2 as thermal fluctuations increases (as can be seen in the bottom panel of Fig. 1 at  T < 1).
Under the influence of magnetic field, the system behavior changes as seen in Figs Albert network influenced by positive magnetic field (Fig. 2a). The field term in the Hamiltonian is effectively a constant holding the network above the mean of two states at 1 2 . As the magnitude of the magnetic field increases, the network takes longer to reach the asymptotic state. For an anti-ferromagnetically coupled modified Ising model of Barabási-Albert network, it can be observed for that for small magnitudes of the positive magnetic field (h << 1), the asymptotic property of order parameter vanishes as in the case of a ferromagnetically-coupled system [ Fig. 2b]. However, for higher magnitudes of the magnetic field, it can be seen that the field term can trigger activity in the network i.e. switch from M = 0 to M = 1 at 0 < T < 1 and subsequently follow the dynamics of a ferromagnetically-coupled system. A negative magnetic field, on the other hand, inverts the dynamics of a ferromagneticallycoupled modified Ising model instead. As can be seen in Fig. 3a, at −2.5 < h < 0 there is an abrupt drop in the order parameter to 0 and for lower values the network remains inactive (as can be verified from our observations in Figs. 2 and 3). An antiferromagnetically coupled network has order parameter M = 0 at h = 0. Lower values of the magnetic field keep the network in the inactive state. For a positive magnetic field, the network undergoes a relatively smooth (almost abrupt) phase transition to the active state. Owing to this, unlike in a ferromagnetically coupled network, it can be observed that intermediate values of order parameter and M → 1 as h increases, confirming our observations in Figs. 2 and 3. Thus it can be inferred that the modified Ising model of a Barabási-Albert network undergoes phase transition due to the magnetic field as shown in Fig. 4. It can be observed that the transition has a discontinuity in order parameter, and hence this may be a first-order phase transition. Hysteresis loops characterize systems that undergo a first-order phase transition. It implies that the network may show more than one value of order parameter for a given magnetic field, h. The hysteresis loop shows the dependence of the state of the system on its history, and it is this phenomenon that forms memory in a hard disk drive.
The procedure to investigate the existence of hysteresis has been well-established, particularly in the context of magnetic materials. We apply the same method for the modified Ising model of a Barabási-Albert network summarized shortly here. Starting with a high negative magnetic field, h, and a stable configuration of the system, we increase the magnitude of the magnetic field slowly. For some value of h, the local field The modified Ising model of a Barabási-Albert network exhibits hysteresis: ferromagnetically coupled, J = 2 (indicated by dots). Simulation parameters: N = 5 × 10 3 and preferentially attached links, m = 5. The gray curve indicates order parameter as the system is driven forward from h 0 = 10 to h n = −10 and the black curve as the system is driven backward from h 0 = −10 to h n = 10 for a node flips. This causes changes in the effective field of the nodes connected to this node, thereby causing them to flip. Once the flipping in the system has thermalized, the order parameter of the system is measured. Subsequently, the magnetic field is increased slightly, and the process repeated until the order parameter attains a stable state. This way, one can obtain one half of the hysteresis loop (for h from −∞ to ∞). The other half of the hysteresis loop is obtained when the magnetic field, h is decreased (for h from ∞ to −∞).
A typical hysteresis loop takes the form of a sigmoid; however, in the case of a ferromagnetically coupled modified Ising model the loop is almost a rectangle, as can be seen in Fig. 5a. We will analyze these observations and discuss the asymptotic behavior in detail using analytical approaches in Sect. 4.

Mean field approximation
One of the most important analytical tool to study disordered systems is represented by mean-field theories. Mean field theory is frequently used due to its conceptual simplicity, as a useful tool, especially when there is no exact solution for the problem. This approximation is used to reduce an interacting problem to a non-interacting one which is easier to solve.

Theorem 1 The critical magnetic field of the modified Ising model of a Barabási-Albert network scales linearly with coupling constant and preferentially attached links used to construct the network.
Proof Let us consider the modified Ising model of a Barabási-Albert network treated numerically in Sect. 3. Rewriting the Hamiltonian of the ferromagnetically-coupled system with gene-type spins 0, 1, where J i j = J A i j . Since the adjacency matrix A i j is symmetric, the factor 1 2 is included so as not to count any pairs twice. We can write the interactions between neighboring spins in terms of their deviations from the average spin M as, Assuming that the fluctuations around the mean spin is small, the Hamiltonian can be rewritten as, Consider the second term in the right hand side of Eq. 5. This can be written as since A i j = A ji , A is symmetric. Therefore from Eqs. 5 and 6, This is the mean-field Hamiltonian for a chosen realization of the network. So the ensemble average of the Hamiltonian of the system is, For a Barabási-Albert network, where k i is the number of links of the ith node of the network Bianconi (2002) (cf. "Appendix"). From Eqs. 8 and 9, using the relation Hence the modified Ising model of a Barabási-Albert network reduces to a system of non-interacting spins in an effective local field, h eff i = (h + J mk i ). The partition function can be evaluated as, The mean spin, M can be calculated from the partition function using the following relation: From this, evaluating ln Z , Therefore from Eqs. 12 and 13, Therefore the central mean-field equation for ferromagnetically coupled Barabási-Albert network with asymmetric spins takes the implicit form, Similarly for anti-ferromagnetically coupled Barabási-Albert network (J → −J ) the central mean-field equation is, Note that the order parameter depends on the coupling constant, J and node degree, k i . Let us first study the behavior of the system in the absence of magnetic field. The meanfield equation for ferromagnetically coupled Barabási-Albert network with gene-type spins and no external field is, where ± stands for ferromagnetically and anti-ferromagnetically coupling respectively. From Eq. 17 we can investigate the asymptotic behavior for ferromagnetically and anti-ferromagnetically coupled modified Ising model of a network. For a ferromagnetically coupled Barabási-Albert network when T → ∞, β J Mk i → 0, so exp(−β J Mk i ) → 1 ⇒ M → 1 2 . As T → 0, β J Mk i → ∞, so exp(−β J Mk i ) → 0 ⇒ M → 1. These confirm the observations in the top panel in Fig. 1. Similarly for an anti-ferromagnetically coupled modified Ising model of a Barabási-Albert network we can verify the limit cases: as On the other hand, as T → 0, exp(β J Mk i ) → ∞ ⇒ M → 0. These validate the observations in the bottom panel of Fig. 1. In order to compare the results of mean-field approximation with Monte Carlo simulations, we have plotted the results using these two different approaches in Fig. 6 which also shows a slight discrepancy between the numerical simulations and the mean-field ansatz. This arises because we neglect the spin product term in the mean-field approximation with the assumption that the fluctuations around the mean spin is small (Eq. 4).
For T >> 1, using Taylor expansion M can be approximated as, M ≈ 1 2±β J m . We can conclude that, for a fixed large T in ferromagnetically coupled systems, those with larger J and M have larger M and vice versa. This investigation predicts the behavior of the system presented in Fig. 7 and validates Monte Carlo simulations. The situation is reversed for an anti-ferromagnetically coupled system due to the presence of plus sign in the denominator. Eqs. 15 and 16 indicates that at T >> 1, However, in both cases, the asymptotic behavior of the system is preserved, for T → ∞ (or β → 0), M → 1 2 (cf. Eq. 18). In the case where T not tending to ∞, the value of M depends on the magnitude and direction of magnetic field, h. This implies that for an anti-ferromagnetically coupled system, when h > J m then m > 1 2 ; however, for h < J m we have m > 1 2 . Similar conclusions can be made when h is negative in a ferromagnetically coupled system. Therefore, the behavior of the system changes at |h c | = J m.
Getting back to our biological model, with J > 0, when the majority of the genes are active above the critical h(h c ), this represents a healthy state. On the other hand, when the majority of cells are inactive below the critical h(h c ), this represents a disease state. This different behaviour is also illustrated in Fig. 2.
In other words, consider the limit of small T where we can neglect fluctuations. In this case, our biological system, have an order parameter 1 for h c > J M and 0 for h c < J M. Therefore, above critical parameter h c , all the genes are active while they suddenly become inactive below h c as a result of a first-order phase transition. This critical external energy interaction depends on the strength of internal interaction (coupling constant J in gene-gene interactions) as well as the parameter m of the Barabási-Albert network.
This approximates the critical magnetic field, h c ≈ 5 for the choice of simulation parameters, which is very close to our observations from numerical simulations as can be verified in Figs. 2 and 3. Although the analytical results predict that the network size does not influence phase transition in the modified Ising model of the Barabási- Albert network, the numerical results predict a weak dependence of h c on network size (Fig. 8c), which appears in systems with large network sizes. The dependence on parameters J and m is over-estimated by the mean-field calculations as can be seen in Fig. 7. In the next Sect. 4.1 we will derive the expression for the critical magnetic This can be mapped to the Hamiltonian of the classical spin system by introducing new spin variables as, For s i = 0 → s i = −1 and for s i = 1 → s i = 1. Substituting the spin variables in the Hamiltonian Eq. 19 we make the H 0,1 → H −1,1 transformation, Since J i j = J ji , N i, j (s i + s j ) = 2 N i, j s i , eq. 21 can be re-written as, So the problem of an Ising model with gene-type spin system is mapped on to a problem of Ising model with classical spin system as, where The fact that even in the absence of magnetic field there is an intrinsic local magnetic field, a N i, j J i j 2 in the system reflects the asymmetricity of the spins present in the problem. In principle, any physical quantity of the system of modified Ising model of a interaction can therefore be derived from the system of Ising spins, However we are interested in the critical magnetic field as derived in Sect. 4.1. Note that, the first term of the right hand of Eq. 23 is a constant and by redefinition of the zero of energy we have, This is the Hamiltonian for a chosen realization of the network. So the ensemble average of the system Hamiltonian is, where, Remark 1 From the above transformation to the classical Ising model, the average critical field for the modified Ising model can be verified.
The average critical field for the system h c can be derived by, where k i is approximated by the average number of links,k. Note thatk = This validates our results presented in Sect. 4.1. Eq. 29 predicts that the critical magnetic field depends linearly on J and m. The numerical simulations confirms the analytical predictions on critical magnetic field ( Fig. 8a and b).

Conclusions
In living systems, collective flipping of coherently expressed genes is associated with disease progression. This flipping causes the step by a step-change in the phenotype of the cell, causing it to transition from normal phase to diseased phase. Similarly, in magnetic systems, collective flipping of spins is associated with the loss of spontaneous magnetization. Therefore it is intuitive to consider gene networks as two-state thermodynamic systems in a heat bath obeying Boltzmann statistics.
On the one hand, in the statistical physics community, there have been extensive studies on phase transitions occurring in Ising models of scale-free networks with classical spins (discussed in Sect. 1). These methods are analytically tractable and applicable to very large sizes. On the other hand, in the systems biology community, a wealth of literature exists that motivates the modeling of networks with binary states for small to medium sizes (discussed in Sect. 2). This work on modified Ising model lies at the intersection of statistical physics and network biology and is presented with an intention to bring the two schools of thought together.
In this regard, we have proposed here an adaptation of a well-established model in statistical mechanics that could be used to study phase transitions in living systems. This model allows simplification of interactions in complex systems; can be studied analytically; and renders itself adaptable to the representation of complex genetic systems, thereby allowing testing of the diverse hypothesis that may cause complex disorders. The modified Ising model is an adaptation of the classical Ising model constructed for networks with a scale-free-like structure and whose activity is described by a binary random variable. This is a general statistical method to deal with poorly understood non-linear large scale models arising in the context of biological networks.
We have presented a basic numerical and theoretical framework to investigate phase transitions using modified Ising model where the control parameters energy and entropy are modeled by the magnetic field and temperature, respectively. Taking the Barabási-Albert model as the toy model, we have shown that such a system undergoes phase transition owing to the influence of the critical magnetic field. This is synonymous to a simple Mendelian disease where there is a strong field near the diseased gene. In complex diseases, the influence of the magnetic field is spread among different genes with different strengths.
The critical magnetic field of the system scales linearly as a function of the number of preferentially attached links and coupling constant. Further, we have shown that the modified Ising model can be mapped to a classical Ising model of a Barabási-Albert network. The simulation setup presented herein can be directly used for any biological network connectivity dataset and is also applicable to other networks that exhibit similar states of activity. The model can be adapted for directed or weighted networks and could also take a continuum of activity states such as in a Potts model. Additional interaction terms may be added to the Hamiltonian to model epistatic interactions between genes. Further, the modified Ising model is capable of predicting the existence of structurally or functionally organized clusters in the network.
We have shown that a purely qualitative model such as the modified Ising model is capable of predicting phase transitions given only the connectivity without logical rules or kinetic parameter data. This indicates that dynamics on these networks may depend more on structure than on the specific details of the processes. The modified Ising model is capable of scaling to networks of sizes up to tens of thousands and can potentially predict similarities between apparently unrelated complex systems.
Funding Open Access funding provided by Projekt DEAL. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through GSC111; and Exploratory Research Space (ERS) Seed Fund 2017 in Computational Life Sciences (CLS001). All simulations were performed using the RWTH Compute Cluster under general use category; priority category allocated to AICES and JRC users; and with specific computing resources granted by RWTH Aachen University under project rwth0348. The authors gratefully acknowledge the generous support of the aforementioned funding and computing resources.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. node j is linked to node i with connectivity k i (t j ) is given by, is proportional to the number of links k i at time t j , and number of preferentially attached links m. The dynamic solution of connectivity at time t i is, From Eqs. 30 and 31 we have, If N is large we can approximate the total number of edges in the network at time t j , given by the sum j α=1 k α as, j α=1 k α = m 0 + 2mt j ≈ 2mt j because m 0 << N . The factor 2 comes from the fact that as we create a link which connects two nodes, the number of links of each of them increases by 1. Substituting Eq. 33 in 32, The adjacency elements of the network A i j are equal to 1 if there is a link between node i and j and 0 otherwise. Consequently the mean over many copies of a Barabási-Albert network From Eq. 31 we can re-write for t = N steps, and similarly, From Eqs. 36 and 37, The average of the adjacency matrix over many realizations can be approximated by the network parameters as,