Rewiring driven evolution of quenched frustrated signed network

A framework for studying the behavior of a classically frustrated signed network in the process of random rewiring is developed. We describe jump probabilities for change in frustration and formulate a theoretical estimate in terms of the master equation. Stationary thermodynamic distribution and moments are derived from the master equation and compared to numerical simulations. Furthermore, an exact solution of the probability distribution is provided through suitable mapping of rewiring dynamic to birth and death processes with quadratic asymptotically symmetric transition rates.


I. INTRODUCTION
Statistical Physics of complex networks has been applied to numerous problems in different areas of science, such as systems biology [1,2], neurobiology [3], sociology [4], economics [5,6], and history [7], just to mention a few.
In the framework of the statistical physics approaches to modeling complex real systems, a particular place belongs to spin-like models, that were naturally extended from the regular lattices to networks.To provide some examples, Ising models have been used to analyze social networks [8] and their tendency to self-organize into communities [9].Potts' models were also used for the same purpose [10,11] as well as in intrusion detection in secured networks [12] and modeling multi-opinion propagation [13] Among all of them, the study of frustrated spin models [14,15] in networks plays a crucial role in understanding and analyzing complex socioeconomic systems.Frustrated models, describing interactions between competing entities [16,17] capture the inherent conflicts, constraints, and interactions that are prevalent in social and economic contexts.Application of frustrated models provides insights into the emergence of collective phenomena [18], decision-making processes, and the dynamics of social networks.These models help elucidate the causes and consequences of socioeconomic imbalances [19], the formation of coalitions and alliances [16], the impact of policy interventions [20], or evolution of welfare [21], to name the few.
Moreover, the study of frustrated models enables the identification of critical points and phase transitions in socio-economic systems, offering valuable information for designing resilient and sustainable solutions to societal challenges.Ultimately, understanding frustrated models allows us to navigate the complexities of socioeconomic systems, predict their behavior, and develop strategies for promoting stability, fairness, and overall societal well-being.Examples of this application include investigations of stable and meta-stable phasesynchronization patterns in neural networks [22], evaluating policies to obtain wanted socioeconomic outcome [20], defining punctuating equilibrium in financial time series [23] etc.Furthermore, their relevance increases when frustrated systems are combined with signed graphs, in which each edge has a positive or negative sign [24].Frustrated models of signed networks have found intriguing applications in deciphering intricate interactions within biological systems.In these systems, nodes represent entities such as genes, proteins, or even organisms, while positive and negative edges signify activation and inhibition relationships respectively.This nuanced approach allows researchers to model the delicate balance of molecular and cellular processes more accurately, unveiling the complexities of signal transduction [25,26], gene regulation [27,28], etc.By integrating both positive and negative interactions, signed networks provide a holistic view of dynamic biological systems, enabling the exploration of feedback loops, bistability, and emergent behaviors.These networks have proven particularly useful in uncovering regulatory mechanisms [27], disease progression [29], and drug targets [30], shedding light on the underlying principles governing life's fundamental processes and offering valuable insights for biomedical research and therapeutic development.
A substantial body of literature delves into the dynamics of signed networks [31][32][33][34], primarily focusing on the dynamics of spins (signs) within these networks, with limited attention given to the evolution of network properties driven by dynamic changes.The pertinence of rewiring algorithms for signed networks finds its most illustrative demonstration in tools like BiRewire3 [35], employed to evaluate the significance of signed network structures in biological research.Nevertheless, analytical outcomes for such algorithms are still outstanding, and the results hinge on intricate simulations of the systems.
Besides the application of frustrated models to realworld problems, their relation to different network parameters is still not solved satisfactorily.The problem, similar to the one analyzed in this paper, was previously addressed in [44].In it, the authors provide an approximate solution to the problem of evolution of ground state energy of frustrated classical Ising system with disorder modeled through the rewiring parameter of the SW model [45].In that paper, the authors give an approximate relation between the energy of frustration and rewiring probability p of the small-world model, based on heuristic arguments.The validity of the relation was supported by simulations, except for large values of parameter p.
In this paper, we comprehensively analyze how the quenched system's count of frustrated links evolves with rewiring, employing an approximate master equation.Furthermore, we demonstrate that the probability distribution of the number of frustrated bonds remains consistent with our calculated values in the thermodynamic limit.

II. ANALYTICAL RESULTS
A frustrated classical Ising system is described in the form of a graph with N nodes, with each node representing a spin that can be in one of two states (+1 or -1).In the terminology used in graph theory, this means that a set of nodes N can be partitioned into two sub-sets of different colors, black (B) and white (W), where each node of spin down, n i = −1 ∈ B is colored black and spin up, n j = 1 ∈ W is colored white.Anti-ferromagnetic interaction between two spins is represented by l-th link, enumerating to |L|.Depending on the state of two spins, we define the sign function σ : l → {−1, 0, 1} which tells us if the chosen link exhibits positive frustration f + i.e. link connects two nodes with positive signs, negative frustration f − i.e. a link connects two nodes with negative signs or is unfrustrated otherwise.As such, we are considering undirected signed networks G = (N, L, σ).
Since our use case is that all the interaction is antiferromagnetic, we can define frustration count F (G) as the number of frustrated links of the graph G depending on the sign function σ: where We start with the network (lattice) that has f + 0 frustrations among positive spins and f − 0 frustrations among negative spins.For the evolution of the system, random rewiring is considered, specifically its most popular variant -the Maslov-Sneppen algorithm [46] which provides a microcanonical evolution of the "smallworld" network.The rewiring of two links is described in table I.We chose this type of rewiring because it preserves degree distribution and can rightfully be considered as a microcanonical process, as opposed to rewiring of only one link per time step, which preserves average degrees but not exact local degrees and can therefore be considered as equivalent to a canonical process.
The main variable we track in the process of rewiring the network is frustration.Only a discrete number of frustration state variable f (t) = f + (t)+f − (t) is created or destroyed in a given moment, and the simplest case is when frustration number changes by the single amount per unit time Here f + (t) and f − (t) represent the number of frustrated links between positive and negative spins at time t.The frustration state of the system also depends solely on the knowledge of the most recent condition and in that case, the Markov assumption is valid, and we can formulate the evolution process in terms of the Master equation.

A. Master equation for rewiring
For integer state space differential Chapman-Kolmogorov equation [47] translates into the Master equation of the form: Initial configuration Ending configuration Change in Frustration Approximate Probability  displaying all potential network rewirings along with their associated probabilities.The initial probability constitutes a basic approximation, as it does not account for the neighboring links.
at the previous time t ′ < t is determined through summation of all other possible k states and corresponding transitional probabilities W (f |f ′ ; t), which can be written as: In our case, using the Maslov-Sneppen rewiring algorithm, in every single time-step the number of frustrated links can increase or decrease only by 2. Therefore we have only two possible transition probabilities, which leads to the general master equation being: Since Maslov-Sneppen rewiring can increase/decrease the number of frustrated links only by 2 with one of them being between positive nodes and the other between negative nodes, this means that the difference between positive and negative frustrated links is conserved in this type of rewiring i.e.
Therefore trivially we have that both f + and f − can be written as We can now specify functions t ± (f ) which describe the probability of change in frustration by consulting the table I.
In the following, we will mainly use the probabilities of choices that are changing the number of frustrated links since they are needed to specify the transition probabilities, while other choices.The probability that frustration decreases means that we choose 1 out of all positive links, and 1 out of all negative links, divided by the number of all combinations of pairs of links.The difference in positive and negative frustration ∆ is a key parameter that distinguishes an ensemble of graphs.
The probability of an increase in frustration via dual rewiring is the number of all unfrustrated pairs divided by the number of all pairs.These approximated probabilities were confirmed by running simulations of choosing two links, rewiring them, and comparing the increase or decrease in frustration.Such an approximation is good enough to compute the evolution of frustrated links subject to the Maslov-Sneppen rewiring.In table I, we provide additional more precise approximations that might be of use for extremely small network systems, but that we did not use in the rest of this study.

B. Stationary solution
If the limit of a large network, the solution of ( 5) is going to converge to simulations.We can integrate this solution directly by numerical integration, using a choice of boundary probabilities p(f < ∆; L, ∆, t) = 0 and p(f > L; L, ∆, t) = 0 with initial condition p(f ; L, ∆, t = 0) = δ(f 0 ).We can also use thermodynamical limit (N → ∞) to approximate the difference between probability distributions at time t + 1 and t as: The usual properties of generating function hold: First, we calculate the simpler case using stationarity condition lim t→∞ ṗ(f, t) = 0 to obtain the solution in the form of a hypergeometric function.
From this expression, we can obtain the probability distribution using an inverse Z-transform.
For f even and in the interval f ∈ [0, L], we obtain that the stationary probability is equal to while vanishes in all the other cases.The expected value of this stationary distribution for L > 0 is simply: When this discrete stationary probability distribution p(f ; L, ∆) is evaluated for large enough parameter L and with ∆ = 0, it tends to normal distribution.

C. Evolution of moments
Though determining the precise distribution p(f, t) is intricate, calculating the expected value f (t) from the differential equation in the generating function is straightforward.The general solution takes the form: same rewiring process for triangular lattice starting from f > 0. The inset axis shows the mean absolute error decreasing for larger system size L.
The value of frustration depends on the parameters f 0 , L, ∆.The stationary state in the limit of a large network is described as ⟨f ⟩ ∼ L 2 + ∆ 2 1 L .In models in which ∆ ∼ O(L) the deviation from the trivial solution will be governed solely by the difference between positively and negatively frustrated bonds.
Although the parameter ∆ influences the final expected value of the evolution, it does not influence the rate with which the system achieves its thermodynamic state.However, for ∆ = 0, we have two initial boundary values: f min and f max , both approaching the same expected value (11) as t → ∞ (as shown in figure 2).

D. Probability distribution evolution
A stochastic birth and death process can be described using the Chapman-Kolmogorov equations: ) where λ n and the µ n are respectively the birth and death rates.In the analysis of our rewiring process, we can use the same equation by defining Since the master equation which governs the process of rewiring always adds or subtracts 2 frustrated links, we introduce variable φ, such that 2φ = f .This leads to the new equation (15) In which Noticing that the smallest number of frustrated links possible in the system is equal to ∆, we re-parametrize it so that δ = ∆/2, and introduce variable ϕ = φ − δ, which leads to parameters: This particular set of parameters corresponds to the birth and date process previously analyzed in the work of Roehner and Valent [48] which solved the general solution for the transition rates: Our process then has parameters defined as: This redefinition of transition rates corresponds to [48] analysis of finite processes where the population can only evolve between the levels 0 and N = L/2.If we take a large limit so that L 2 ≫ L The general formula derived is: where [a] = the integer part of a and H l being: Details of the computation are provided in the appendix.

III. VALIDATION
To test the analytical solutions, numerical simulations were executed exploiting the Python Networkx package.
Our main objective is to validate that our estimated analytical model aligns with the outcomes of simulations for two scenarios: the transition from a state with 1) minimal frustration and 2) maximal frustration.In the second scenario, it's important to note that this isn't the stationary state of maximum conceivable frustration, but rather a state where a change in frustration value is plausible, i.e. the most likely outcome.

A. Minimum and maximum Frustration
First, we randomly assign spin values to the nodes of a graph of size L to determine the ∆ parameter.We then optimize the graph to minimize frustration and obtain the system configuration.Starting from the configuration set by parameters (L, f 0 , ∆ 0 ), we proceed with the rewiring process for t = 3000 steps.At each step t, we calculate the system's frustration, repeating this process for n = 1000 iterations.The mean value and standard deviation are obtained by averaging n frustration values: In the case of a square grid, achieving a minimum frustration value of f = 0 is always possible.However, for a triangular grid, the minimum frustration is not zero.We can observe that the theoretical solution approximates the rewiring process effectively.Accuracy improves with larger system sizes, evident from the reduced mean absolute error between the analytical function and the average frustration value.The systems' frustration tends towards the expected theoretical value.
Clearly, a system that has f = L, ∆ = L won't exhibit the behavior of maximum possible decrease in frustration by process of rewiring noted in figure 2.
For a system to be in the state of maximum frustration and allow for this type of evolution, two conditions need to be satisfied, namely that f 0 = L and that the difference between positive and negative frustrated links ∆ = 0. To achieve that, a simulated annealing algorithm was used with a Euclidean distance cost function that targets maximum frustration and minimum ∆.
Where ŷ = (L, 0) is the target value of the system state.This condition would ideally mean that the graph has two separate parts with f + + f − = L.Not all graph configurations can satisfy these two conditions, and for them, the next optimum solution is chosen, an example is shown in figure 4.
The value of average frustration decreases as predicted in figure 5.The ratio f /L □ = 0.89, 0.85, 0.81 and f /L △ = 0.90, 0.95, 0.86 for square lattice and triangular lattice respectively denote the initial starting points.Because of the stochastic behavior of simulated annealing and the implication of the size of the system, in practice, these ratios are slightly lower for larger systems.To calculate the probability distribution of a system, we have several options, one of which is to assume that the process is ergodic after a certain period of rewiring, arriving at the stationary state.We count the values of frustration for each of the n systems, therefore obtaining the probability distribution, and take the average over the stationary period ∆t = t 2 − t 1 .This is shown on the right plots of figures 9 and 10 (in appendix) for the square and triangular lattice of similar size for increase and decrease of frustration by the process of rewiring, respectively.
Red lines in figures 9 and 10 (in appendix) denote the means and the standard deviation calculated from the analytic distribution which matches the simulated values nicely, suggesting that, although the theoretical FIG.5: Rewiring process starting from the approximately maximum obtainable frustration and lowest ∆ value for square and triangular lattice.
model is only approximate, the extra information contained in the probability distribution describes the process well for larger systems.
A second way to compare the simulated probabilities with analytical results, depending on the time in the rewiring process, is to sample all n evolutions in k groups and count the occurring number of all possible frustrations (f = 0..., L) in each of the k systems for every time step t.In simulations we mainly used k = 100.We can then calculate the average probability distribution ⟨p(f, t)⟩ and its standard deviation.
To measure the similarity between the theoretical stationary distribution and the simulated one, we use Jensen-Shannon divergence: FIG.6: Jensen -Shannon Divergence between asymptotic analytical p(f ) distribution and simulated average distribution q(f, t).
Where D(p||q) is the Kullback-Leibler divergence: In figure 6 we see that this metric measure decreases from the value of 1 representing the difference between the certain probability of system in state f 0 and the stationary distribution to a smaller nonzero value, implying that the distributions become similar as the rewiring proceeds.
In order to demonstrate that our analytical approach can follow the evolution of the whole probability distribution, we present figure 7. The initial steep increase of the Jensen-Shannon divergence is related to the discrete nature of time steps for the simulations, which gets smoothed out during the evolution as the number of time steps becomes larger.
Here an important Remarque is in order.Assuming that the stationary probability p(f ; L, ∆) is the probability that a given system has frustration f , these values do not depend on the time t of the rewiring process.If we sample n = 1000 systems with random configurations of ∆ we can also obtain the probability distribution of frustration for a system of size L. It is clear that the stationary probability distribution obtained by the process of rewiring will coincide with this probability distribution, as is shown in figure 8. FIG.7: Jensen-Shannon divergence of the simulated and analytically evaluated probability distribution at time t.Initial jump is related to the initial difference between evolutions mainly governed by the difference between discrete time steps of simulated probability distributions and the continuous time steps of the master equation solution.
FIG. 8: Probability of a system assuming frustration value f in the stationary regime depending on the size of the system L. Simulation is represented in color, while the analytic solution is black.

IV. CONCLUSION
In this paper, we presented an estimated solution for the energy of a frustrated Ising antiferromagnetic model based on the disorder level in the system.Though earlier studies offered some approximate solutions, our work introduces a more accurate solution.Additionally, we present supporting evidence that demonstrates the precise agreement between our solution and an exact solution in the thermodynamic limit.
From a purely physical perspective, we have investigated the quenched spins, in a fast-changing network as opposed to another limit of slowly changing network and fast flipping spins studied in [44].While the quenched frustration energy changes with ∼ e − (2L−1) L(L−1) t , the previ-ously studied ground state energy per spin changes with t.An interpolating behavior between these 2 regimes that were not previously studied is of particular interest.We envision analyzing such a cross-over regime in the future.Computed evolution of the frustration subjected to random rewiring provides other researchers with important information that can be used in a number of different applications.For example, provided solutions could be used in testing algorithms in extremely large networks.Furthermore, although we focused our attention on frustrated links in signed networks, a very similar analysis can be used in networks with differently colored edges, like for example in multiplex networks [34,[49][50][51][52].We believe that it can be especially useful in the application of frustrated models in large changing net-work systems, such as social, biological, or technological networks.Red full and slashed lines represent the mean value and standard deviation calculated from the theoretical solution.

FIG. 1 :
FIG. 1: Example of the signed lattice graph with described coloring.Negative frustrated links are colored red and positive frustrated links are colored green depending on the node colors.

FIG. 2 :
FIG. 2: Evolution of mean value of frustration through the process of rewiring from boundary values of initial frustration f min = 0 and f max = L for L = 84, 180, 420, with ∆ = 0.

FIG. 3 :
FIG. 3: (a) Rewiring process for a square lattice starting from the lowest possible frustration f = 0 (b) same rewiring process for triangular lattice starting from f > 0. The inset axis shows the mean absolute error decreasing for larger system size L.

FIG. 4 :
FIG. 4: (a) A plot of (f, ∆) configurations.Blue points represent the randomly initialized configuration sample of 1000 graphs, red trajectory shows the simulated annealing process from starting point (84, 84) to an optimum solution (b) Solution to optimization of a L = 84 triangular graph with f = 76, ∆ = 0

FIG. 10 :
FIG. 10: (a) Rewiring process of square lattice for f 0 = 0, L = 760, ∆ = 0 (b) Probability distribution function averaged over time period t 1 = 2000, t 2 = 3000 with theoretical stationary probability distribution shown in black.Red full and slashed lines represent the mean value and standard deviation calculated from the theoretical solution.

TABLE I :
Table