2.1 Dynamical systems modeling and structural balance theory in international relations

Dynamical systems models have been useful in understanding and predicting the dynamics of complex systems in a diversity of fields such as biology, engineering, and sociology (Strogatz & Dichter, Reference Strogatz and Dichter2016). While statistical models are useful in classifying and predicting phenomena, statistical models often fail to provide explicit mechanistic explanations (Midlarsky, Reference Midlarsky1984; Thompson, Reference Thompson1986; Roberts et al., Reference Roberts, Augustine, Lawton, Lindsay, Thiele, Izquierdo and Lockery2016); dynamical systems models, on the other hand, can elucidate the mechanisms responsible for the observed underlying phenomena (Trotta et al., Reference Trotta, Bullinger and Sepulchre2012; Morrison & Kutz, Reference Morrison and Kutz2021). Dynamical systems models of the interaction between states or armed sub-state actors such as insurgents can help identify mechanisms for system destabilization as well as mechanisms for reinforcing stability. These models can be used to determine the risk for minor, localized conflicts to escalate and spread extensively to other nodes in the network. Moreover, they provide a rich set of diagnostic tools for the analysis of escalating conflicts which result in the outbreak of widespread wars.

The application of dynamical systems modeling to international relations has largely been built upon Richardson’s pioneering model of arms races (Gleditsch, Reference Gleditsch2020; Stauffer, Reference Stauffer2021). However, this stream of research predated the surge of scholarship on political networks and so did not engage with questions of network structure, although the Richardson model itself is amenable to network formulation (Ward, Reference Ward2020). Nor did it leverage the concepts and tools of modern bifurcation theory.

The application of network analysis to international relations has bloomed over the past two decades. Most of this work has involved unsigned networks in which either cooperative or, less frequently, conflictual ties are treated as separate networks between nations or militant groups (Hafner-Burton et al., Reference Hafner-Burton, Kahler and Montgomery2009; Maoz, Reference Maoz2011; Victor et al., Reference Victor, Montgomery and Lubell2016; Dorussen et al., Reference Dorussen, Gartzke and Westerwinter2016; Zech & Gabbay, Reference Zech and Gabbay2016). Signed networks, however, allow for an integration of cooperation and conflict consistent with the fact that alliances are usually formed with an eye toward a potential foe and decisions to militarily confront an opponent typically consider who might come to its aid. The investigation of world politics using signed networks has overwhelmingly been through the lens of structural balance theory (Cartwright & Harary, Reference Cartwright and Harary1956), a strand of research that precedes by decades the recent network turn in international relations study.

In its basic form, structural balance theory is formulated on (undirected) triads and asserts that there are only two balanced, and hence enduring, triads: one with all positive ties (“the friend of my friend is my friend”) and one with a single positive tie and two negative ones (“the enemy of my enemy is my friend”). The two other possibilities—two positive ties and one negative or three negative ties—are imbalanced and hence unstable. A network as a whole is considered balanced if it is comprised of only balanced triads. These rules imply that a complete network can have only two perfectly balanced states: a harmonious state in which all ties are positive and the two hostile factions state in which the ties within a faction are positive while the ties between factions are negative (Cartwright & Harary, Reference Cartwright and Harary1956).

There has been a long-running back-and-forth about whether the international system is characterized by structural balance (Harary, Reference Harary1961; Healy & Stein, Reference Healy and Stein1973; McDonald & Rosecrance, Reference McDonald and Rosecrance1985; Maoz et al., Reference Maoz, Terris, Kuperman and Talmud2007) with more recent work tending to support an overall tendency toward balance (Lerner, Reference Lerner2016; Kirkley et al., Reference Kirkley, Cantwell and Newman2019; Burghardt & Maoz, Reference Burghardt and Maoz2020). Going beyond the question of whether the international system displays a tendency toward balance, its level of balance has been observed to fluctuate considerably over time (Doreian & Mrvar, Reference Doreian and Mrvar2019; Burghardt & Maoz, Reference Burghardt and Maoz2020), which suggests that other forces beyond structural balance may be at work. Evidence for structural balance processes has also been found in friendship networks (Kirkley et al., Reference Kirkley, Cantwell and Newman2019), wild mammals (Ilany et al., Reference Ilany, Barocas, Koren, Kam and Geffen2013), and gang networks where inter-gang violence increases among gangs in imbalanced triads (Nakamura et al., Reference Nakamura, Tita and Krackhardt2020).

2.2 Modeling structural balance dynamics

Given a network of initial ties, dynamical systems models can be formulated that evolve ties over time consistent with structural balance theory. The simplest formulation of structural balance dynamics evolves a symmetric network of signed edge weights, $\textbf{X}(t) \in \mathbb{R}^{N \times N}$ , as the following system of coupled, nonlinear differential equations:

(1) \begin{align} \frac{dX_{ij}}{dt} = \sum _{k=1}^N X_{ik}X_{kj} \end{align}

If nodes $i$ and $j$ have either both positive ties or both negative ties with a third node $k$ , the product $X_{ik}X_{kj}$ will be positive and act so as to increase their tie value $X_{ij}$ . In contrast, if their ties with $k$ are of opposite sign, the effect will be to decrease $X_{ij}$ . Both of these effects act to increase the balance in the system. The sum determines whether the net effect of all mutual ties increases or decreases $X_{ij}$ .

Although this equation appears to have been first proposed by Lee et al. (Reference Lee, Muncaster and Zinnes1994), it was not investigated in that paper. It was proposed independently by Kulakowski et al. (Reference Kulakowski, Gawronski and Gronek2005) who found that in simulations this model leads to perfectly balanced outcomes consistent with the static theory. Marvel et al. (Reference Marvel, Kleinberg, Kleinberg and Strogatz2011) used the matrix equation formulation of Equation (1),

(2) \begin{align} \frac{d\textbf{X}}{dt} = \textbf{X}^2 \end{align}

to obtain an analytical solution for $\textbf{X}(t)$ . They showed that the network will generically evolve into a balanced end state as determined by the eigenvector of the initial adjacency matrix $\textbf{X}(0)$ having the most positive eigenvalue, which grows the fastest. If this first eigenvector consists of all positive components, then the system will evolve into the harmonious state. On the other hand, if the first eigenvector has components of both positive and negative sign, then the system will evolve into the two hostile factions state. Marvel et al. (Reference Marvel, Kleinberg, Kleinberg and Strogatz2011) used this model to show good agreement with the alliances which formed in World War II and the split of the Zachary karate club. Morrison & Gabbay (Reference Morrison and Gabbay2020) extended this model to initial conditions containing community structure generated by a two-block stochastic model. They found that phase transitions in the initial network structure determine whether the two final factions align with the block structure or instead have random memberships unrelated to the blocks or, alternatively, whether it is the harmonious state that emerges.

Equation (1) can be analyzed for the case of a directed network in which case the system evolves to four factions instead of two from random initial conditions (Veerman, Reference Veerman2018). It is also possible to define a variant model for directed networks that replaces $X_{ik}X_{kj}$ on the right-hand side of (1) with $X_{ik}X_{jk}$ , which can be shown to evolve into two factions for random initial conditions (Traag et al., Reference Traag, Dooren and Leenheer2013). Although our focus here is on deterministic dynamical systems, structural balance theory can also be implemented as a stochastic model (Antal et al., Reference Antal, Krapivsky and Redner2006).

As we employ a control theory framework to investigate the patterns of network perturbations that are most destabilizing, we note some previous work applying control theory to structural balance. Gao & Wang, (Reference Gao and Wang2018) have investigated the origins of balancing forces with respect to nodal interactions and have considered how to control structural balance edge dynamics using node dynamics, with the assumption that the dynamics of nodes and edges are interdependent. Wongkaew et al. (Reference Wongkaew, Caponigro, Kułakowski and Borzì2015) have considered how to control structural balance dynamics by introducing a “leader” to the system, while Summers & Shames (Reference Summers and Shames2013) have evaluated the control abilities of existing nodes in a network by changing ties.

While the simple model of Equation (1) has yielded great insight into structural balance dynamics and will serve as the main point of departure, it also has important shortcomings. Although it is able to determine the outcomes that emerge from a network solely under structural balance dynamics, it has no equilibrium, evolving ties toward positive or negative infinity (and does so in finite time (Marvel et al., Reference Marvel, Kleinberg, Kleinberg and Strogatz2011)). Once the model is turned on, the network perforce evolves into a state of complete war or harmony and so it cannot capture the stable, pre-escalatory state which characterizes the international system for long stretches of time before and after systemic wars. Its empirical application therefore depends upon the assumption that structural balance dynamics are either off or on and so, while useful for predicting the composition of the opposing sides once the war starts, the model is silent about what patterns of relationships might set off such a war in the first place. Nor can the model describe de-escalation dynamics, from systemic war back to a peaceful state since it inherently seeks to increase balance. Accordingly, we aim to build upon the dynamical systems model of conflict escalation in Equation (1) to create a model that can manifest a wider range of behaviors in the international system—stability in a peaceful, low-conflict state, transitions to systemic war, and de-escalation.

4.1 Bifurcations in simulation

Figure 1(a) and (b) compare the tie dynamics for a network governed by pure structural balance dynamics (Equation 1) versus a network governed by the simultaneous operation of the direct dyadic and structural balance forces (Equation (3) with $X_{Tij}=0$ ). As noted above, the pure structural balance dynamics results in ties that blow up to positive or negative infinity. This phenomenon occurs for all initial conditions with no means of reverting back to a peaceful state with low tie values.

Figure 1. Simulations of signed network model dynamics. (a) Time series of tie values $X_{ij}(t)$ for pure structural balance model, Equation (1), showing unbound dynamics. (b) Top panel. Time series for the dyadic force and structural balance model, Equation (3), as the balance sensitivity $\alpha$ is changed at the times indicated by the dotted lines; $\beta = 1$ , $L=8$ , $N=10$ , $X_{Tij}=0$ , and $X_{ij}(0)=X_{Dij}$ . $X_{Dij} \sim 0.8 \mathcal{N}\,(0,1) +0.4 N u_{Ci} u_{Cj}$ is a random symmetric matrix with polarized community structure, where $\textbf{u}_C$ is the vector that generates the block structure (see Section 6). Middle panel. Snapshots of dynamic network tie matrices at $t=10$ , $30$ , and $50$ . Note that the color scales are different for each matrix. Bottom panel. The standard deviation of network ties, $\sigma (\textbf{X})$ , over time. (c) Standard deviation of network ties as a function of the structural balance parameter $\alpha$ as well as the initial state (war or peace). The vertical dotted lines are the predicted critical values of $\alpha$ for the peace-to-war and war-to-peace bifurcations from Equations (6) and (20) respectively.

In Figure 1(b), the dyadic bias parameters, which also serve as the $t =0$ tie values, are randomly set but have an underlying two-block community structure. Roughly speaking, in signed networks communities are characterized by a greater density of positive ties within communities than between them and, conversely, a greater density of negative ties between communities than within them (Traag et al., Reference Traag, Doreian and Mrvar2019). In two-block structure, the nodes are ordered so that each of the two blocks consists of contiguous nodes. If the ties within each block clearly tend to be more positive than the ties between blocks, then the blocks correspond to distinct communities. Note that it is possible for the blocks to comprise separate communities even when both intra and inter-block ties tend to be positive as long as the former are denser or more positive. If, furthermore, the intra-block ties tend to be positive while the inter-block ones tend to be negative then the community structure of the bias network is said to exhibit “polarization” in that there is a hostile relationship between the two communities.Footnote ⁴ Both the bias and dynamic networks can exhibit polarization, but the dynamic network, through its evolution by the model, has the potential to greatly exacerbate the polarization built into the bias network (the bias parameters themselves are fixed during the model simulation).

Figure 1(b) shows the dynamics of Equation (3) in which the structural balance sensitivity, $\alpha = L/\gamma$ , is modified at times $t=15$ , $30$ , and $45$ (by changing $\gamma$ ). The initial state is given by the bias network, which is taken to have polarized structure. At the parameter shift times, however, the dynamic network state is not reset to $\textbf{X}_D$ , but rather the present tie values serve as the initial conditions for the new $\alpha$ value and the system is allowed time to reach equilibrium.

The first time segment from $t=0$ to $t=15$ with $\alpha =0.05$ shows that the system quickly reaches an equilibrium in which the tie values are close to their initial $X_{Dij}$ ones. Thus, unlike pure structural balance dynamics, our model possesses a stable state spanning a range of relatively small, positive and negative tie values as seen in the low standard deviation of the ties in the bottom plot of Figure 1(b), which corresponds to the peace state defined in the preceding section. It accords with recent work that seeks to provide a finer-grained characterization of peace as shades of rivalry and friendship, rather than as an undifferentiated state of simply “not war” as has often been assumed in conflict studies (Diehl, Reference Diehl2016). Observe that the polarization of the dynamic network in the peace state (and hence the bias network) is reflected in the $\textbf{X}(10)$ matrix plot by the yellowish-green color of the diagonal blocks signifying overall positive intra-community ties and the bluish color of the off-diagonal blocks signifying overall negative inter-community ties.

In the second time segment from $t=15$ to $t=30$ , $\alpha$ is doubled to $0.1$ and the peace state destabilizes and the system transitions to the war state. The ties diverge into two bunches of strong positive and negative tie values, yielding a much larger standard deviation than the peace state. The increase in structural balance has greatly increased the polarization in the dynamic network as seen in the $\textbf{X}(30)$ matrix plot: not only are the intra and inter-block ties uniformly positive and negative respectively, but their magnitudes are much higher as well. At $t=30$ , $\alpha$ is lowered back down to $0.05$ , its value in the first time segment. Yet, this does not bring the system back to the peace state. $\alpha$ must be decreased even further, to $0.01$ at $t=45$ , to induce the transition back to the low-deviation peace state. This behavior is an example of hysteresis, where a system’s future state is determined not only by its parameter values but also by its current state (Strogatz & Dichter, Reference Strogatz and Dichter2016).

Figure 1(c) shows the tie standard deviation, $\sigma (\textbf{X})$ , for various $\alpha$ values and initial conditions. For $\alpha \lt 0.018$ , only the peace state is stable and the dynamic network will converge to it regardless of whether the initial condition corresponds to the peace or war state. For $\alpha \in (0.018, 0.053)$ , both the peace and war states are stable and the system is bistable as is consistent with hysteresis; which of the two equilibria the network converges to is dependent on initial conditions. For $\alpha \gt 0.053$ , only the war state is stable, and the network will converge to it for both peace and war state initial conditions.

4.2 Bifurcations in the eigenvalue spectrum

It will be helpful to analyze the matrix of the dynamic network and its time evolution in terms of its eigenvalues and eigenvectors. The set of eigenvectors $\textbf{s}_i$ of the real, symmetric matrix $\textbf{X}$ satisfy the following equation:

(4) \begin{align} \textbf{X} \textbf{s}_i = \lambda _i \textbf{s}_i \end{align}

where the eigenvalues $\lambda _i$ are real and ranked in order of descending value. The eigenvectors $\textbf{s}_i$ form an orthonormal set when the eigenvalues are distinct and $\textbf{X}$ is a full-rank matrix. We make these assumptions in our analysis. In Section 5.1, however, consider a special case in which $\textbf{X}$ is not full rank. In matrix form, the eigenvector decomposition can be written as

(5) \begin{align} \textbf{X} = \textbf{S} \Lambda \textbf{S}^T, \end{align}

where the columns of $\textbf{S}$ are the $\textbf{s}_i$ and $\Lambda$ is a diagonal matrix containing the corresponding eigenvalues $\lambda _i$ . The eigenvectors form an alternative, and often more intuitive, coordinate basis for $\textbf{X}$ . In unsigned networks, for example, the first eigenvector components are proportional to the node eigenvector centralities whereas, for signed networks, strong two-faction structure will be reflected in the first eigenvector in which the factional memberships are of opposite sign.

Figure 2 compares the simulation results in which the matrix of the bias network, $\textbf{X}_D$ , is either completely random with no structure at all or contains polarized community structure, where as noted previously, polarization is a subset of community structure marked by two blocks where positive ties form preferentially within blocks and negative ties between them. Polarized structure can be generated and studied using stochastic block models (Section 6). The dyadic biases also serve as the initial conditions $\textbf{X}(0)=\textbf{X}_D$ for each simulation run. In Figure 2(a), we plot the adjacency matrix spectrum as $\alpha$ is increased. Each horizontal slice contains all the eigenvalues of $\textbf{X}$ at equilibrium for a given $\alpha$ value. The peace-to-war bifurcation occurs at the critical value of the $\alpha$ parameter, $\alpha ^*_{P \rightarrow W}$ . For $\alpha \lt \alpha ^*_{P \rightarrow W}$ , the system is stable in the peace state and the eigenvalues form a dense band although the first two eigenvalues show some small growth as $\alpha$ is increased from zero. For $\alpha \gt \alpha ^*_{P \rightarrow W}$ , the peace state is destabilized and the system transitions to the war state. The leading eigenvalue is now much larger than the others, making the equilibrium network close to rank 1. The curves in Figure 2(a-b) show that the analytically computed leading eigenvalue in the peace and war states (Equations (17) and (19)), and the bifurcation value for $\alpha$ (Equation (6)) are in good agreement with the simulation results. These approximations are good despite the small value for $L$ , which makes such approximations less robust. Figure 2(c) visualizes the equilibrium dynamic network at three $\alpha$ values as well as the corresponding time series of the tie weights. At low $\alpha$ values, the network ties remain close to the random initial conditions set by $\textbf{X}_D$ , while beyond $\alpha ^*_{P \rightarrow W}$ the final network is polarized into two camps with large positive ties within each group and large negative ties between the two groups (the nodes could be reordered so the communities appear as contiguous blocks).

Figure 2. Effect of increasing structural balance sensitivity, $\alpha$ on eigenvalues, $\lambda _i$ , and balance levels, $\eta$ , of the dynamic network at equilibrium. $L=2$ , $\beta =1$ , and $N=10$ . (a) Stable state eigenvalues as a function of $\alpha$ resulting from a random bias network containing no community structure, $\textbf{X}_{Dij} \sim 0.8 \mathcal{N}\,(0,1)$ . Theoretically computed bifurcation value $\alpha _{P \rightarrow W}^*$ (Equation (6)), $\lambda _P$ (Equation (17)), and $\lambda _W$ (Equation (19)). (b) Stable state eigenvalues resulting from a bias network with polarized community structure, $X_{Dij} \sim 0.8 \mathcal{N}\,(0,1) +0.8 N u_{Ci}u_{Cj}$ . (c) Equilibrium dynamic networks and tie time series for $\alpha = 0.1$ , $\alpha = 0.07$ , and $\alpha = 0$ for $\textbf{X}_D$ without community structure. (diagonals set to zero) (d) Equilibrium dynamic networks and tie time series for $\alpha =0.1$ , $\alpha =0.03$ , and $\alpha =0$ for $\textbf{X}_D$ with community structure. (e) Balance levels $\eta$ of the equilibrium network (blue circles) as a function of $\alpha$ for the no initial community structure case. Lower and upper range of $\eta$ values from null model simulations shown as gray lines. (f) Equilibrium network balance levels for the initial community structure case. The green curves in (a) and (b) show the analytical expression, (17), for the first eigenvalue in the peace state.

The spectrum of the case where $\textbf{X}_D$ contains polarized community structure is shown in Figure 2(b). For $\alpha =0$ , there is no structural balance force and the tie values remain at their dyadic biases as seen in the time series of Figure 2(d), so the equilibrium $\textbf{X}$ spectrum is the same as $\textbf{X}_D$ . Unlike the no structure case, the first eigenvalue is substantially larger than the rest reflecting the underlying two-block structure of $\textbf{X}_D$ , which can be seen in the $\alpha =0$ matrix visualization. This community structure results in a much smaller $\alpha ^*_{P \rightarrow W}$ than the no structure case. This suggests that bias networks exhibiting clear polarized structure, strong enough to appear in the first eigenvalue, may be more easily destabilized than those without. We present a more general argument about this behavior in Section 6.

Structurally balanced networks are comprised of balanced triads—those containing either one or three positive ties. Balance levels in signed networks can be measured by counting the number of balanced versus imbalanced triads that exist in the network. Kirkley et al. (Reference Kirkley, Cantwell and Newman2019) develop a global balance metric $\eta$ , where $\eta = 0$ corresponds to perfect balance and $\eta = 1$ corresponds to the average level of balance generated in null model simulations where the null models are generated by randomly swapping the signs of the ties in the observed network. Figure 2(e–f) shows the balance level of the network at its final equilibrium state relative to the range of balance levels observed in the null model instances, shown by the gray bands. When the blue dot moves outside of the gray band this indicates that the system is significantly more balanced than expected polarized the null process. The no initial community structure case does not show a significant level of balance until after the bifurcation (Figure 2(e)) at which point the dynamic network transitions to the war state and becomes perfectly balanced. The case where the bias network has initial community structure, in contrast, shows significant levels of balance, which smoothly increases with $\alpha$ even before the bifurcation (Figure 2(f)). These plots show that the onset of war greatly increases the level of structural balance in the dynamic network. Yet, the community structure plot also reveals that it is possible for the dynamic network to be balanced in a statistical sense in the peace state. The fact that structural balance tends to characterize peace and increases due to major war has been observed empirically for the international system as noted above (See Section 2).

The peace-to-war bifurcation threshold is given by

(6) \begin{align} \alpha ^*_{P\rightarrow W} = \frac{\beta }{4 \lambda _{D1}}, \end{align}

where $\lambda _{D1}$ is the leading eigenvalue of $\textbf{X}_D$ . We derive this expression in Appendix A. The values for $\alpha ^*_{P\rightarrow W}$ predicted by this formula are in good agreement with the observed thresholds as shown by the red dashed lines in Figure 2(a) and (b). The community structure case has a higher $\lambda _{D1}$ than the no structure case as seen in the $\alpha =0$ spectra, which by (6) lowers the bifurcation threshold.

Even without the derivation, Equation (6) can be motivated intuitively. That increasing the direct dyadic force strength $\beta$ increases $\alpha ^*_{P\rightarrow W}$ makes sense as the dyadic force seeks to maintain the peace state of low conflict and cooperation. When the elements of $\textbf{X}_D$ are drawn from a zero mean distribution as in Figure 2, the squared eigenvalue, $\lambda _{D_i}^2$ , is proportional to the variance of $\textbf{X}_D$ carried by each eigenvector.Footnote ⁵ The appearance of $\lambda _{D1}$ in the denominator of (6) then reflects the fact that networks with greater values of $\lambda _{D1}$ are more polarized in the peace state and so more unstable to war.

We note that for sufficient net positivity in $\textbf{X}_D$ , a bifurcation between the peace state and the harmony state of universal high cooperation occurs as the balance sensitivity $\alpha$ is increased, analogous to the peace-to-war transition. It is also possible to have a direct transition between the war and harmony states. In the pure structural balance dynamics model, Equation (1), the war-to-harmony transition occurs when strong two-faction community structure becomes subordinate to harmonious one-faction structure in the initial network (Morrison & Gabbay, Reference Morrison and Gabbay2020). We consequently expect our model (3) to likewise undergo this bifurcation for sufficiently strong structural balance and bias network structure. To illustrate this bifurcation, instead of varying $\alpha$ , we will vary the dyadic bias parameters via the outgroup affinity $p_{\text{out}}^+$ , the (conditional) probability of forming a positive tie with the outgroup in the stochastic block model which can be used to generate $\textbf{X}_D$ (see Section 6). In Figure 3(a-b), the reversal of dominant structure is observed in the crossing in the spectrum of $\textbf{X}_D$ and in the qualitative change in its leading eigenvector from polarized structure to a homogeneous structure. Recall that $\textbf{X}_D$ plays a double role as the initial value of the dynamic network and as persistent parameters in the dyadic bias force. For low $p_{\text{out}}^+$ values the leading eigenvector contains two factions while for large $p_{\text{out}}^+$ values the leading eigenvector contains homogeneous structure (Figure 3(b)). The spectrum of the final network, however, does not display a crossing (Figure 3(c)); nor is there a discontinuity in the first eigenvalue as in Figure 2. Instead, it is the components of the first eigenvector which manifest the discontinuous signature of the bifurcation (Figure 3(d)). The block structure of the first eigenvector components below $p_{\text{out}}^+=0.5$ yields the war state (Figure 3(b)), whereas the uniformly positive components above $p_{\text{out}}^+=0.5$ yields the harmony state.

Figure 3. Simulation illustrating war-to-harmony bifurcation. (a) Eigenvalue spectrum of $\textbf{X}_D$ , also equivalent to initial dynamic network $\textbf{X}(0)$ , for varying outgroup affinity $p_{out}^+$ . Sample matrices at left depict how the ties are more positive and community structure is weaker above the crossing at $p_{out}^+=0.5$ than below. (b) Components of leading eigenvector of $\textbf{X}_D$ . When the leading eigenvector only has values $\geq 0$ , the corresponding leading eigenvalue is labeled as “homogeneous” while when the leading eigenvector has both positive and negative values, the corresponding eigenvalue is labeled as “factional”. (c) Eigenvalues of the equilibrium network after simulation of Equation (3) with $N=100$ , $\beta = 1$ , $\alpha = 0.02$ , and $L=4$ . (d) Components of leading eigenvector of equilibrium network. Other stochastic block model parameters (see Section 6): $d_{in}=d_{out} = 0.4$ , $p_{in}^+=1$ .

5.1 Special case: Equal and opposite factions

Consider the special case where there are two equally sized, opposing factions in which all dyadic bias parameters have the same magnitude $\mu$ , but are positive for ties within the same faction and negative between factions so that $X_{Dij}=\pm \mu$ . We set the initial condition to be $X_{ij}(0)=X_{Dij}$ . Equation (3) then reduces to two equations,

(7) \begin{align} \frac{dx}{dt} = -\beta (x-\mu )+L \tanh\!\left(\frac{\alpha N}{L}x^2\right) \end{align}

(8) \begin{align} \frac{dy}{dt} = -\beta (y+\mu )-L \tanh\!\left(\frac{\alpha N}{L}y^2\right) \end{align}

for ingroup ties $x$ and outgroup ties $y$ , with initial conditions, $x(0)=\mu$ and $y(0)=-\mu$ . However, if we substitute $y=-x$ in the bottom equation, we recover the top one and its initial condition, showing that the two equations are not independent. The system therefore reduces to one equation as solving for $x$ using (7) then determines $y$ ; the conflictual inter-faction ties are simply the negative of the cooperative intra-faction ties.

We make Equation (7) more tractable by replacing the hyperbolic tangent term on the positive axis with a piecewise linear approximation consisting of a line of unit slope connecting the origin to a plateau of height $L$ such that $L^{\ast} \tanh(\theta/L) \rightarrow \theta$ for $0 \leq \theta \lt L$ and $L^{\ast} \tanh(\theta/L) \rightarrow L$ for $\theta \geq L$ , where $\theta =\alpha Nx^2$ . This approximation for the structural balance term is shown by the purple dashed line in Figure 4 (note it is piecewise linear in $\theta$ not $x$ , which is why the transition section is concave upwards). Under this approximation, Equation (7) becomes

(9) \begin{align} \frac{dx}{dt} = -\beta (x-\mu )+ \alpha Nx^2, \qquad 0 \leq x \lt \sqrt{\frac{L}{\alpha N}} \end{align}

(10) \begin{align} \frac{dx}{dt} = -\beta (x-\mu )+ L,\qquad x \geq \sqrt{\frac{L}{\alpha N}} \end{align}

Using these equations, we can derive simple formulas for the fixed points in the system: the peace and unstable fixed points will be obtained from (9) and the war fixed point from (10).

Figure 4. Strength of structural balance (purple) and direct dyadic (blue) forces and behavior of the equilibria in the special case. The approximations of the structural balance force (dashed) used in (9) and (10) agree well with the exact form (7) (solid), but are less accurate in the shoulder region. Stable (solid circles) and unstable (open circles) fixed points occur at intersections of balance and dyadic forces. $\alpha =0.03$ , $L=6$ , $N=10$ , and $\mu = 1/2$ .

The equilibria of the above equations are found by setting $dx/dt=0$ , which implies that the direct dyadic and structural balance forces are of equal magnitude but of opposite sign. The right-hand side of Equation (9), being quadratic, has two roots. One yields the peace equilibrium, $x_p$ , which is stable and occurs at

(11) \begin{align} x_P = \frac{\beta - \sqrt{\beta (\beta - 4 N \alpha \mu )}}{2 N \alpha } \end{align}

The other fixed point is unstable and occurs at

(12) \begin{align} x_U = \frac{\beta + \sqrt{\beta (\beta - 4 N \alpha \mu )}}{2 N \alpha } \end{align}

Setting the right-hand side of Equation (10) to zero yields the fixed point corresponding to the war state, $x_w$ which is stable:

(13) \begin{align} x_W = \frac{L}{\beta } + \mu \end{align}

The above fixed points refer to the ingroup tie values—the outgroup tie values are their negatives. We observe that the structural balance force enters into these fixed points differently: for the war state via its maximum value, $L$ , whereas for the peace and unstable states via the sensitivity $\alpha$ , which affects how much the structural balance force grows in response to changes around its neutral point. Equation (13) also explains our identification of the war tie value scale with $L/\beta$ (see Section 3). Since the peace tie value scale is $\mu$ , our assumption that the war tie values are much higher than the peace ones is equivalent to $L/\beta \gg \mu$ .

Figure 4 shows how the nature of the equilibrium solutions changes as the dyadic force strength parameter $\beta$ is varied while the structural balance parameters remain constant. The fixed points are located at the intersections of the linear dyadic force and the nonlinear structural balance force, where their magnitudes are equal but point in opposite directions with the dyadic force acting to decrease $x$ (toward $\mu$ ) and the balance force acting to increase $x$ . The stability of a fixed point with respect to small perturbations is determined by the relative magnitudes of the forces on either side of it. At moderate levels of the dyadic strength, $\beta = 1$ , there are three fixed points in the system. The war equilibrium at $x_W\approx 6$ is stable since perturbations to its left (right) will be pulled back by the stronger balance (dyadic) force. The peace equilibrium at $x_P \approx 0.6$ is similarly stable. But the intermediate fixed point, $x_U$ is unstable since the stronger balance force pulls positive perturbations away from $x_U$ and towards $x_W$ and negative perturbations are drawn toward $x_P$ by the stronger dyadic force so that $x \rightarrow x_W$ , when $x\gt x_U$ and $x \rightarrow x_P$ , when $x\lt x_U$ . When the dyadic strength is strong, $\beta = 2$ , the stable peace state in the only fixed point, so that $x \rightarrow x_P$ for all $x$ . Finally, when the dyadic stability is weak, $\beta = 1/2$ , the stable war state at $x_W \approx 12.5$ is the only fixed point so that $x \rightarrow x_W$ for all $x$ (although it is difficult to discern in the figure, the dyadic and balance forces do not intersect at low $x$ ).

The bifurcation responsible for the changes between the different solution regimes is a saddle-node bifurcation. Two fixed points disappear when a stable fixed point collides with the unstable one as a parameter is changed. For example, if we increase $\beta$ starting from the $\beta =1$ condition of Figure 4, the intersections determining $x_U$ and $x_W$ will draw closer together until they coincide when the dyadic force is tangent to the shoulder of the exact structural balance force (or intersects the corner of the approximate one). This coincidence point marks the war-to-peace bifurcation as it disappears upon further increase in $\beta$ , leaving only the $x_P$ fixed point. Analogously, the peace-to-war bifurcation occurs when $x_P$ and $x_U$ collide, leaving only $x_W$ .

We can now derive expressions relating to the parameters at the bifurcation points. The peace-to-war bifurcation occurs when $x_P = x_U$ . Equating (11) and (12) implies that the square root term must vanish, which yields the condition (ignoring the $\beta =0$ solution, which is simply pure structural balance),

(14) \begin{align} \alpha ^*_{P\rightarrow W} = \frac{\beta }{4 N \mu } \end{align}

where we have solved for the critical value of the structural balance sensitivity. Similarly, the war-to-peace bifurcation occurs when $x_U=x_W$ . Equating (12) and (13) yields

(15) \begin{align} \alpha ^*_{W\rightarrow P} = \frac{L}{N(\mu +L/\beta )^2} \end{align}

Equation (14) shows that the peace-to-war transition occurs at larger $\alpha$ values for larger dyadic force strengths, $\beta$ , a behavior that accords with the dyadic force’s role in stabilizing the peace equilibrium. Increasing the dyadic bias parameter, $\mu$ , destabilizes the system at a lower $\alpha ^*_{P\rightarrow W}$ , signifying that greater polarization in the peace state makes it more susceptible to war. Increasing the number of nodes $N$ also causes the peace state to become more fragile.Footnote ⁶ The dependence on $L$ on the right-hand side of Equation (15) implies that the war-to-peace bifurcation sees the height of the balance force plateau, unlike the peace-to-war transition. Since $\alpha ^*_{W\rightarrow P}\rightarrow 0$ as $L\rightarrow \infty$ , a larger $L$ , which corresponds to a larger level of cooperation and conflict in the war state, makes it more difficult for the system to transition back to the peace state.

5.2 General results

We consider bifurcations in the dynamics for the general case where the dyadic bias parameters can be different for each dyad. In addition, we allow for application of a transient perturbation which can induce a transition from the peace to the war state (or more generally a state of high tie values). We state the key results in this section and show the derivations leading to them in Appendix A.

Our dynamical system of edges, Equation (3), can be written succinctly as a matrix dynamical system

(16) \begin{align} \frac{d\textbf{X}}{dt} = -\beta (\textbf{X} - \textbf{X}_D) + L \tanh \!\left ( \frac{\alpha \textbf{X}^2}{L} \right ) + \textbf{X}_T(t) \end{align}

where $\textbf{X}$ is the symmetric matrix of ties $X_{ij}$ , $\textbf{X}_D$ is the matrix of the dyadic bias parameters, $X_{Dij}$ , and $\textbf{X}_T(t)$ is the matrix of the transient perturbations $X_{Tij}(t)$ . We can analyze the stability of the dynamics in the low-dimensional eigenvalue-eigenvector space of $\textbf{X}$ and derive the system’s bifurcation conditions in terms of the leading eigenvalue of $\textbf{X}_D$ .

We first state the results in the absence of an imposed perturbation $\textbf{X}_T(t)$ . As in the special case, the peace-to-war transition is a saddle-node bifurcation. It occurs when the respective eigenvalues, $\lambda _P$ and $\lambda _U$ , of the peace and unstable equilibria, $\textbf{X}_P$ and $\textbf{X}_U$ , merge. The general analysis approximates the structural balance term using $\tanh \theta \approx \theta$ for small $\theta$ as in the special case. The leading eigenvalue of the peace state is given by

(17) \begin{align} \lambda _P = \frac{\beta - \sqrt{\beta ^2 - 4\alpha \beta \lambda _{D1}}}{2 \alpha } \end{align}

where $\lambda _{D1}$ is the leading eigenvalue of the dyadic bias matrix $\textbf{X}_D$ . The lowest eigenvalue of the unstable state is the same as the expression for $\lambda _P$ , but with a positive sign in front of the square root term (see Equation (A18)). Setting $\lambda _P=\lambda _U$ , equivalent to the vanishing of the square root term above, yields the critical value for the structural balance sensitivity in the war-to-peace bifurcation,

(18) \begin{align} \alpha ^*_{P\rightarrow W} = \frac{\beta }{4 \lambda _{D1}} \end{align}

After the peace-to-war bifurcation, the system converges to the war equilibrium, $\textbf{X}_W$ , which has leading eigenvalue $\lambda _W$ . Using the flat plateau approximation of the structural balance force, $L \tanh \theta \approx \pm L$ for large $|\theta |$ yields the following approximation for $\lambda _W$ ,

(19) \begin{align} \lambda _W = \lambda _{D1} + \frac{L N}{\beta } \end{align}

The bifurcation from the war state back to the peace state occurs via a different saddle-node bifurcation where $\lambda _U = \lambda _W$ . The approximate condition for the war-to-peace bifurcation is then given by

(20) \begin{align} \alpha ^*_{W\rightarrow P} &= \frac{NL}{(\lambda _{D1} + NL/\beta )^2} \end{align}

The general expressions for the critical $\alpha$ values, (18) and (20), are the same as in the special case except that $N\mu$ has been replaced by $\lambda _{D1}$ (the next section shows that $\lambda _{D1} = N\mu$ in the special case). Figures 1 and 2 demonstrate how the analytically computed bifurcation values match the bifurcations observed through simulation. However, the peace-to-war estimated bifurcation is more accurate than the war-to-peace estimated bifurcation.

The analysis of the system with the transient perturbation term is restricted to the case when $\textbf{X}_T(t)$ is applied as an impulse in which all of its elements are constant for the duration of the impulse. In this case, the impact of $\textbf{X}_T(t)$ can be accounted for by a modified dyadic bias matrix, $\tilde{\textbf{X}}_D$ , that combines $\textbf{X}_D$ and $\textbf{X}_T$ via

(21) \begin{align} \tilde{\textbf{X}}_D(t) = \textbf{X}_D + \frac{\textbf{X}_T(t)}{\beta } \end{align}

The threshold condition for $\textbf{X}_T(t)$ to induce the system to switch from the peace to the war state is analogous to the peace-to-war bifurcation in which $\lambda _P=\lambda _U$ . However, the leading eigenvalue of the dyadic bias matrix in Equation (17) for $\lambda _P$ as well as that for $\lambda _U$ (Equation (A18)) is replaced by the leading eigenvalue of $\tilde{\textbf{X}}_D$ , which we write as $\tilde{\lambda }_{D1}$ . Accordingly, the threshold condition is the same as the peace-to-war bifurcation condition, Equation (18), but with $\tilde{\lambda }_{D1}$ instead of $\lambda _{D1}$ . For the imposed perturbation context, however, we keep the model parameters fixed and vary $\textbf{X}_T$ , so it is more appropriate to cast the condition in terms of the critical value of $\tilde{\lambda }_{D1}$ ,

(22) \begin{align} \tilde{\lambda }_{D1}^* &= \frac{\beta }{4 \alpha } \end{align}

7.1 Minimum energy destabilization using dyadic bias leading eigenvector

The perturbation to the peace state that requires the least amount of energy, $||\textbf{X}_T||$ , to destabilize is along the direction of the leading eigenvector of $\textbf{X}_D$ , denoted by $\mathbf{s}_{D1}$ . The resulting perturbation matrix can be written as $\textbf{X}_T =\sigma \mathbf{s}_{D1} \mathbf{s}_{D1}^T$ . We solve for the minimum value of $\sigma$ that will destabilize the system ( $||\textbf{X}_T||=\sigma$ ). Destabilization requires that the leading eigenvalue of the modified dyadic bias matrix, $\tilde{\lambda }_{D1}$ , is at least as large as the critical eigenvalue so that

(24) \begin{align} \tilde{\lambda }^*_{D1} & \leq \tilde{\lambda }_{D1} \end{align}

Making use of Equation (22) and the definition of $\tilde{\textbf{X}}_{D}$ , (21), then gives

(25) \begin{align} \frac{\beta }{4 \alpha } \leq \boldsymbol{\lambda _1} \!\left (\textbf{X}_D + \frac{\textbf{X}_T}{\beta } \right ) \end{align}

(26) \begin{align} \leq \boldsymbol{\lambda _1} \!\left (\textbf{X}_D + \frac{\sigma }{\beta } \mathbf{s}_{D1} \mathbf{s}_{D1}^T \right ) \end{align}

(27) \begin{align} \leq \lambda _{D1} + \frac{\sigma }{\beta } \end{align}

where the notation $\boldsymbol{\lambda _1}(\!\cdot\!)$ refers to the first eigenvalue of its argument and we have used $\textbf{X}_D \textbf{s}_{D1} = \lambda _{D1} \textbf{s}_{D1}$ and $\mathbf{s}_{D1}^T \mathbf{s}_{D1} = 1$ to show that $\textbf{s}_{D1}$ is also an eigenvector of $\tilde{\textbf{X}}_D$ . Solving for $\sigma$ we find the condition required for destabilization,

(28) \begin{align} \sigma \geq \beta \!\left(\frac{\beta }{4 \alpha } - \lambda _{D1}\right) \end{align}

It is also useful to write this condition as

(29) \begin{align} \sigma \geq \beta (\tilde{\lambda }^*_{D1} - \lambda _{D1}) \end{align}

7.2 Destabilizing directions for perturbations

Although $\textbf{X}_T = \sigma \textbf{s}_{D1} \textbf{s}_{D1}^T$ is the perturbation that destabilizes the system with the minimum amount of energy (minimum $||\textbf{X}_T||$ ), there are other, slightly higher energy perturbations that can also destabilize the system. We can search for perturbation directions within the space of perturbations $||\textbf{X}_T||=\sigma$ that come closest to destabilizing the system. These are the perturbation directions that are most likely to destabilize the system since they require less energy to destabilize the system than other perturbations. We use the following optimization scheme to find perturbation directions $\textbf{X}_T$ that are closest to destabilizing the system for a given energy level $||\textbf{X}_T|| = \sigma$ :

(30) \begin{align} \text{Energy minimizing scheme:} \quad \textbf{X}_T^{op}=&\mathop{\text{argmin}}_{\textbf{X}_T \in A}\left [ ||\tilde{\lambda }_{D1}^* - \tilde{\lambda }_{D1}|| \right ] \end{align}

where $A = \{\textbf{X}\;:\; ||\textbf{X}|| = \sigma \}$ , $\tilde{\lambda }_{D1}$ is the leading eigenvalue of $\tilde{\textbf{X}}_D$ , and $\tilde{\lambda }_{D1}^*$ is the stability threshold (Equation (22)). Energy minimizing perturbations result in minimal changes, on average, in tie strengths, whether positive or negative, since $||\textbf{X}_T||$ has been minimized. Equation (30) brings the system close to destabilization by perturbing the system in a direction close to the leading mode. We describe the destabilization scheme in more detail in Appendix B and Figure B1.

We now present two alternative destabilization schemes. The low-rank war state structure that appears after destabilizing the peace state will not necessarily contain two equally sized factions. The leading eigenvector could instead contain all nodes or mostly all nodes of the same sign. The interpretation of destabilization in this case would be instead the formation of a tightly-knit alliance structure where, due to structural balance dynamics, nodes are compelled to form positive ties with other nodes due to mutual allies. We can use the “eigenvector polarization” measure of the leading eigenvector, $\phi _1$ , to determine the extent to which the leading eigenvector contains community structure or a lack of community structure close to homogeneity (Morrison & Gabbay, in preparation). The eigenvector polarization $\phi _i$ of eigenvector $\textbf{s}_i$ of the adjacency matrix $\textbf{X}$ is $\phi _i = \textbf{s}_i^T \textbf{M} \textbf{s}_i$ , where $\textbf{M}$ is the signed modularity matrix (Traag et al., Reference Traag, Doreian and Mrvar2019; Newman, Reference Newman2006). The modularity matrix clusters nodes into communities.Footnote ⁸ If $\phi _1$ is the largest eigenvector polarization, then the network structure is dominated by polarization between two hostile factions.

We use $\phi _1$ to define the following two schemes, which seek to: (1) bring the system close to destabilization with a widespread two-faction conflict (Equation (31)), or, alternatively, (2) bring the system close to destabilization with a widespread alliance bloc (Equation (32)). $\textbf{X}_T^F$ indicates a perturbation that promotes factional structure in the leading eigenvector. A destabilizing perturbation that promotes factional structure will induce polarization in the system. $\textbf{X}_T^H$ indicates a perturbation that promotes homogeneous structure in the leading eigenvector. A destabilizing perturbation that promotes homogeneous structure will induce harmony in the system.

(31) \begin{align} \text{Polarizing:} \quad \textbf{X}_T^F=&\mathop{\text{argmin}}_{\textbf{X}_T \in A}\left [ ||\tilde{\lambda }_{D1}^* - \tilde{\lambda }_{D1}|| + \varepsilon \frac{1}{\phi _1} \right ] \end{align}

(32) \begin{align} \text{Harmonizing:} \quad \textbf{X}_T^H=&\mathop{\text{argmin}}_{\textbf{X}_T \in A}\left [ ||\tilde{\lambda }_{D1}^* - \tilde{\lambda }_{D1}|| + \varepsilon \phi _1 \right ] \end{align}

Here, $\phi _1$ is the first eigenvector polarization measure of $\tilde{\textbf{X}}_D$ , which changes as the perturbation $\textbf{X}_T$ is varied, and $\varepsilon$ is its weighting. Forcing large changes in ties, whether increasing positive relations or increasing hostility requires more “energy” (a larger $||\textbf{X}_T||$ ) than smaller perturbations in relations between countries. Because the leading eigenvector of the adjacency matrix may not contain factions, energy minimizing perturbations may not result in factions. The polarizing perturbation includes a faction-promoting term $\frac{1}{\phi _1}$ , weighted by $\varepsilon$ , that is minimized when $\phi _1$ is large. This promotes a solution where the leading eigenvector of the system contains factions. The extent to which factions are promoted is determined by the strength of $\varepsilon$ (small enough values for $\varepsilon$ may not produce solutions that contain factions). The harmonizing perturbation, in contrast, deincentivizes factions in the leading eigenvector by including the faction-suppressing term $\phi _1$ . The homogeneous perturbation is minimized when $\phi _1$ is small which promotes a solution where the leading eigenvector of the system does not contain factions. These perturbation schemes assume that $\textbf{X}_T$ can be made large enough to result in a transition from the peace state to an equilibrium with high conflict and cooperation tie values (Figure B1(b)). If the maximum structural balance force strength, $L$ , is sufficiently large, then the final state of the system will be determined by the first eigenvector of $\tilde{\textbf{X}}_D$ (see Equation (A24)). Therefore, as $\textbf{X}_T$ can cause the leading eigenvector of $\tilde{\textbf{X}}_D$ to be very different from that of $\textbf{X}_D$ , the high tie value state produced by the perturbation may contain very different factions than those in the war state arising from the peace-to-war bifurcation. While Equation (31) promotes factions in the systemic war equilibrium state and Equation (32) promotes a large alliance bloc encompassing most, if not all of the nodes, Equation (30) does not favor a particular structure in the war state, it simply emphasizes the structure most present in the dyadic bias matrix.

8.1 World War I onset due to triggering perturbations

We use the maximally destabilizing directions formulas, Equations (30), (31), and (32), to find destabilizing directions for perturbations $\textbf{X}_T$ for the network of great powers leading up to the First World War under the dynamics of Equation (16). World War I is considered a classic example of a systemic war that emerged from the complex dynamics of five European great powers, the United Kingdom, France, Germany, Austria-Hungary, and Russia. We will compare our computationally derived destabilizing perturbation directions to ties that have been historically identified as critical to the outbreak of the Great War.

In order to illustrate our method on a simple system, we chose to only include the five most powerful European countries in the year 1913 in our network, as measured by the “Composite Index of National Capability” (CINC) score (Singer et al., Reference Singer, Bremer and Stuckey1972; Singer, Reference Singer1987). The scores range from $0.16$ (Germany) to $0.06$ (Austria-Hungary). We exclude countries below this threshold as we do not weight the ties by country power and including the ties of less relevant countries adds “noise” to the relevant network structure. The next most powerful country below the cutoff is Italy with a score of $0.03$ , half that of Austria-Hungary. The political science literature also identifies these five countries as the most relevant to the outbreak of WWI (Kagan, Reference Kagan1995; Taylor, Reference Taylor1954). In other work we consider larger networks with many more countries (Morrison & Gabbay, in preparation); we are able to incorporate countries of vastly different sizes into the same network in this work by weighting each country’s ties by country power, resulting in a network that accurately reflects the relevance of each tie.

We initialize the system using empirical data as follows. $\textbf{X}_D$ is set to the tie values in the year 1913 using the militarized interstate disputes (MIDs), alliances, and rivalries data from the Correlates of War dataset and the Handbook of International Rivalries (Palmer et al., Reference Palmer, D’Orazio, Kenwick and Lane2015; Gibler, Reference Gibler2009; Singer & Small, Reference Singer and Small1966; Thompson & Dreyer, Reference Thompson and Dreyer2011). Alliances contribute positive weights, MIDs contribute negative weights for states on opposite sides of a dispute and positive weights for states on the same side of a dispute, and rivalries contribute negative weights. Weights between states are added and scaled to range from −2 (rivals) to + 2 (allies), making all values in the dyadic bias matrix low-conflict/cooperation in comparison with the war levels of $\pm 5$ as is roughly consistent with our assumption of peace and war having different tie scales.Footnote ⁹ (Table 2). Self-tie values are set to zero in the bias network. The dyadic force strength is set to $\beta = 1$ , the balance sensitivity to $\alpha = 0.03$ , and $L=5$ . With these parameter values the nonlinear dynamics contain two stable equilibrium, the peace state and war state, with the peace state sitting close to the unstable equilibrium. With initial conditions close to $\textbf{X}_D$ , the network dynamics stabilize close to the bias network tie values. Since $\textbf{X}_D$ contains values in the range $(\!-\!2,2)$ , the network stabilizes at the peace state when no perturbation is added to the system. The size of $\alpha$ does not affect the perturbation direction that is most destabilizing but instead the magnitude of perturbation necessary to surpass the destabilization threshold. We assume the system is sufficiently close to the peace-to-war bifurcation so that it resides in the bistable zone where both the war and peace equilibria are stable and that it can be destabilized with perturbations smaller than the dyadic bias matrix, $||\textbf{X}_T|| \lt ||\textbf{X}_D||$ . Setting $\alpha = 0.03$ accomplishes these objectives. Making $\alpha$ significantly larger would move the threshold past the necessity for a perturbation and the peace state would destabilize without the need for a perturbation. Setting $\alpha$ significantly smaller would indicate that a destabilizing perturbation larger than the dyadic bias matrix would be required to destabilize the system and so perturbations would be unlikely to disrupt the status quo.

The left panel of Figure 5(a) shows the dyadic bias matrix, $\textbf{X}_D$ while the middle panel plots the eigenvalues of $\textbf{X}_D$ , the eigenvalues of the modularity matrix $\textbf{M}$ associated with $\textbf{X}_D$ , and the associated eigenvector polarization $\phi$ of all the eigenvectors. The leading eigenvalue of $\textbf{X}_D$ has value $\lambda _{D1}= 6.01$ , which is observed to be well separated from the remainder of the spectrum as is the first eigenvalue of the modularity matrix. The fact that $\phi _1$ is the largest eigenvector polarization signifies that $\textbf{X}_D$ ’s first eigenvector $\textbf{s}_{D1}$ contains contrast-like community structure consisting of hostile factions (Morrison & Gabbay, in preparation). The two factions present in the leading eigenvector are the United Kingdom (UKG), France (FRN), and Russia (RUS) versus Germany (GMY) and Austria-Hungary (AUH). Also displayed is the critical eigenvalue of the modified dyadic bias matrix for destabilization, calculated from Equation (22) to be $\tilde{\lambda }_{D1}^* = 8.33$ . The optimum perturbation $\textbf{X}_T^{op}$ that enables destabilization with the minimum energy is displayed in the right panel of Figure 5(a). It is proportional to the matrix formed by the first eigenvector of $\textbf{X}_D$ and is given by $\textbf{X}_T^{op} =\sigma \mathbf{s}_{D1} \mathbf{s}_{D1}^T$ as explained in Section 7.1. The outer product of the leading eigenvector shows that FRN, GMY, AUH, and RUS are strongly entrenched in their respective factions while the outer product values for UKG are much closer to zero, indicating a looser affiliation with its faction. The minimum energy is calculated from Equation (29), which yields $\sigma = 2.32$ .

Figure 5. (a) $\textbf{X}_D$ is set to the tie weights of the great powers in 1913, immediately preceding WWI (left). Eigenvalue spectrum of $\textbf{X}_D$ and modularity matrix, and eigenvector polarization spectrum (middle). Note that $\phi_1$ is rightmost eigenvector polarization. For $\alpha = 0.03$ and $\beta =1$ , the bifurcation threshold is $\tilde{\lambda }_{D1}^* = 8.33$ . The leading eigenvalue of $\textbf{X}_D$ is $\lambda _{D1}=6.01$ resulting in $\sigma = 2.32$ . Minimum energy perturbation $\textbf{X}_T^{op}$ (right) leads to final state below. (b) Three perturbation directions requiring minimal energy (Equation (30)). Perturbations applied in these directions result in two factions in the leading eigenvector, and therefore two factions in the final state, $\textbf{X}(t \rightarrow \infty )$ . The factional structure perturbations (Equation (31)) are the same as the minimum energy perturbations for the WWI network. (c) Maximally destabilizing directions that have harmonizing structure (Equation (32)) result in a leading eigenvector with an almost homogeneous structure. Error minimized, $E = \tilde{\lambda }_{D1}^* - \tilde{\lambda }_{D1}$ . Equation (A24) approximates equilibrium war states. Diagonal values in equilibrium states manually set to zero in network visualizations.

Substantively, the optimal perturbation matrix, $\textbf{X}_T^{op}$ , is marked by two tight factions—Germany and Austria-Hungary versus France and Russia; the United Kingdom is loosely aligned with the Franco-Russian bloc. This structure is consistent with the antagonisms and alliances that were most entwined in the several years before the war (Taylor, Reference Taylor1954; Kagan, Reference Kagan1995) Russia and Austria-Hungary competed for influence in the Balkans while France and Germany had severe territorial disputes and a history of animosity. France and Russia were allies as were Germany and Austria-Hungary. It is therefore reasonable that a crisis which activated tensions between the members of these alliances would be most efficient in triggering a war. In contrast, although the United Kingdom was aligned with France and Russia as a member of the Triple Entente, it was reluctant to get involved in Continental conflicts and its participation in the event of a war was uncertain.

The rightmost panel of Figure 5(b) shows the equilibrium state, $\textbf{X}(t\rightarrow \infty )$ , to which our dynamical model evolves the system as a result of a transient perturbation in the direction of $\textbf{X}_T^{op}$ . The war state reached by the system matches the two opposing sides present during WWI with the United Kingdom, France, and Russia in one faction and Germany and Austria-Hungary in the other.

Going beyond the least energy perturbation, $\textbf{X}_T^{op}$ , we solve for alternative, sparse perturbation directions using the optimization functions delineated above, setting $\sigma = 2.32$ , the energy level for the optimal solution. We compute the error as the amount to which the leading eigenvalue deviates from the stability threshold $E = \tilde{\lambda }_{D1}^* - \tilde{\lambda }_{D1}$ . The error tells us how close the system with a perturbation $||\textbf{X}_T||=\sigma$ in the given direction is to destabilizing. Smaller values mean that the system is close to a bifurcation.

Figure 5(b) shows three optimal sparse perturbation directions, $\textbf{X}_T$ , computed using Equation (30) from 100,000 randomly generated perturbations. These sparse perturbations have a substantial component in the direction of $\textbf{X}_T^{op}$ . All three involve only the Continental powers. The middle one, $\textbf{X}_{T2}$ , depicts flare-ups in the two core rival dyads, France vs. Germany and Russia vs. Austria-Hungary. It was the assassination of Austrian Archduke Franz Ferdinand by a Serbian nationalist that set off a crisis between Austria-Hungary and Russia, Serbia’s ally, which ultimately led to the war. The right-hand perturbation, $\textbf{X}_{T3}$ , reflects increased cooperation between Germany and Austria-Hungary and between France and Russia as well as escalating tensions within the Austro-Hungarian–Russian and German–French rivalries. During the crisis, both Germany and France promised to come to the aid of their respective allies in the event of war whereas Britain, absent from this perturbation, was much more ambiguous about its intentions (Kissinger, Reference Kissinger1994). The lowest energy one, $\textbf{X}_{T1}$ , represents an alternative scenario in which the conflicting dyads are France vs. Austria-Hungary and Russia vs. Germany with, again, support between allies. As with the optimal perturbation in Figure 5(a), all three perturbations when evolved by the dynamical model converge to the same two factions as present during WWI, as seen on the rightmost matrix of Figure 5(b).

We use Equation (31) to compute sparse perturbations that promote polarized structure and find the same perturbations as those found using the minimum energy scheme, and hence the polarizing scheme is not displayed as the results are identical to Figure 5(b). This occurs because the leading eigenvector of $\textbf{X}_D$ is already factional. These sparse perturbations have a substantial component in the direction of $\textbf{X}_T^{op}$ . If the leading eigenvector of $\textbf{X}_D$ does not contain factional structure, then the perturbations obtained via Equation (31) will differ from the perturbations obtained via Equation (30). However, since the leading eigenvector of the WWI matrix does contain factional structure, Equations (31) and (30) return similar perturbations. As with the optimal perturbation, the perturbations in Figure 5(b) show that the United Kingdom’s ties with Continental Europe do not appear as significant in triggering conflict. This aligns with historical accounts that the United Kingdom was reluctant to get involved in Continental conflicts (Taylor, Reference Taylor1954; Kagan, Reference Kagan1995).

Figure 5(c) shows the harmonizing perturbations $\textbf{X}_T^H$ computed using Equation (32). The factional structure is much more dominant than other structures in $\textbf{X}_D$ resulting in a larger perturbation required to obtain a harmonious final state, $||\textbf{X}_T^H|| = 4\sigma$ . Ideally, the harmonizing perturbations would result in a homogeneous structure among all nodes in the final equilibrium state. With a low level of perturbation energy, however, the system only manages to positively coalesce four of the five countries in a positive alliance bloc. Perturbations $\textbf{X}_{P1}^H$ and $\textbf{X}_{P3}^H$ lead to a homogeneous coalition against Germany while $\textbf{X}_{P2}^H$ leads to a homogeneous coalition against Austria-Hungary.

Because the network, in this case, contains relatively few nodes ( $N=5$ ), self-ties become a more relevant part of the dynamics. If self-ties are set to zero initially, the self-tie contribution to the force, $\alpha (X_{ii} + X_{jj})X_{ij}$ , is also initially zero, yet grows in relevance as the self-ties grow, reinforcing the sign and strength of each tie. This adds to the escalation and makes it easier for the system to pass the destabilization threshold. We can also consider adding self-ties to the dyadic bias matrix, for example $\textbf{X}_D \rightarrow \textbf{X}_D + c I$ , where $c I$ is the identity matrix scaled by $c$ . In this case the spectra of $\textbf{X}_D$ is shifted by $c$ , moving the system closer to the peace-to-war bifurcation. In both of these cases, self-ties contribute to system instability.

8.2 World War I onset due to changes in dyadic bias community structure

We now take a different approach and evaluate which individual ties between two countries are most critical for increasing the leading eigenvalue and polarization of the dyadic bias matrix, $\textbf{X}_D$ , using tie values from the year 1913, preceding the outbreak of WWI. Given that increasing the first eigenvalue and polarization makes the system more susceptible to war, this analysis can identify which dyads most destabilize the system (or stabilize it if the first eigenvalue decreases). In Figure 6 we systematically shift the ties of $\textbf{X}_D$ by $\varepsilon =\pm 1$ and observe the change each edge perturbation produces in the leading eigenvalue $\lambda _1$ and polarization $\phi _1$ , of $\textbf{X}_D$ . Despite $\textbf{X}_T$ consisting of only a single edge perturbation in each experiment, this small perturbation applied to $\textbf{X}_D$ can result in large shifts in the leading eigenvalue and eigenvector polarization measure (Figure 6(a-b)). For example, Figure 6(b) illustrates how increasing hostilities between France and Austria, $\varepsilon _{2,4} = -1$ , results in an increase in the leading eigenvalue. We represent each edge perturbation’s effect on $\lambda _1$ and $\phi _1$ in a succinct matrix format in Figure 6(c-f).

Figure 6. (a) Great powers network in 1913. (b) Increasing hostilities between France and Austria increases the leading eigenvalue of the network connectivity matrix. (c) Changes in $\phi _1$ for a $\varepsilon = -1$ edge perturbation applied individually to each pair of nodes. Colors for each box in the matrix indicate the change in $\phi _1$ when a perturbation is applied to that pair of nodes. (d) Changes in $\phi _1$ for $\varepsilon = +1$ perturbations. (e) Changes in $\lambda _1$ for $\varepsilon = -1$ perturbations. (f) Changes in $\lambda _1$ for $\varepsilon = +1$ perturbations. Self-ties were not considered.

Figure 6(c-f) shows that the tie between France and Germany as well as the tie between Germany and Austria-Hungary were the most crucial in increasing polarization in the network. The ties between Great Britain and the Continental countries appear to be the least important, implying that while Great Britain formed an important part of the network, changes in ties with this state least impact the stability of the system. Other significant ties include FRN-AUH, FRN-RUS, GMY-RUS, and AUH-RUS. This finding is to be expected as the German alliance with Austria-Hungary and their disputes with Russia as well as the alliance between Russia and France were all factors that proved integral to the development of World War I (Kagan, Reference Kagan1995; Taylor, Reference Taylor1954).