# Steering cooperation: Adversarial attacks on prisoner’s dilemma in complex
networks

Kazuhiro Takemoto  Department of Bioscience and Bioinformatics, Kyushu
Institute of Technology, Iizuka, Fukuoka 820-8502, Japan
<EMAIL_ADDRESS>
###### Abstract

This study examines the application of adversarial attack concepts to control
the evolution of cooperation in the prisoner’s dilemma game in complex
networks. Specifically, it proposes a simple adversarial attack method that
drives players’ strategies towards a target state by adding small
perturbations to social networks. The proposed method is evaluated on both
model and real-world networks. Numerical simulations demonstrate that the
proposed method can effectively promote cooperation with significantly smaller
perturbations compared to other techniques. Additionally, this study shows
that adversarial attacks can also be useful in inhibiting cooperation
(promoting defection). The findings reveal that adversarial attacks on social
networks can be potent tools for both promoting and inhibiting cooperation,
opening new possibilities for controlling cooperative behavior in social
systems while also highlighting potential risks.

  *

Keywords: prisoner’s dilemma, evolutionary game, complex networks, adversarial
attacks

##  1 Introduction

The prisoner’s dilemma game is a fundamental model for investigating the
evolution of cooperation among self-interested individuals [1, 2, 3]. This
game presents a social dilemma where mutual cooperation maximizes collective
benefit, yet defection is each player’s rational choice. This paradigm
reflects universal challenges in various real-world scenarios, such as
resource management, environmental conservation, and collective action, making
the search for solutions highly significant.

Initially, research on the prisoner’s dilemma game in complex networks focused
on lattice networks with uniform connectivity, revealing that spatial
structure promotes cooperation evolution [4, 5, 6]. Recent advancements in
network science [7, 8] have uncovered non-trivial connectivity patterns in
real-world social networks (e.g., small-world topology [9] and heterogeneous
or scale-free connectivity [10, 11]). Subsequent studies explored the impact
of these realistic network structures on cooperation evolution. Heterogeneous
connectivity often promotes cooperation due to hubs [12, 13, 14], although
some studies challenge this view [15, 16]. Conversely, small-world networks
tend to inhibit cooperation emergence due to shortcuts or long-range
interactions [9, 17, 18].

Further investigations examined the impact of link weights on the evolution of
cooperation. Both heterogeneity [19, 20, 21, 22, 23] and adaptive changes in
link weights [24, 25, 26, 27] can promote cooperation. Proposed approaches
include strengthening high-payoff links [24], adjusting weights based on
payoff comparisons [25], and reinforcing links between cooperators [27]. This
concept relates to promoting cooperation through adaptive link rewiring [28,
29, 30].

While these studies provide valuable insights, developing practical methods
for intentional cooperation promotion remains challenging due to the costly
nature of significant network alterations. The concept of adversarial attacks,
originally developed for neural networks [31, 32, 33, 34, 35], may offer a
solution. Recent applications of adversarial attacks to network dynamics,
particularly in voter models, have shown that strategically applied small
perturbations can potentially reverse vote outcomes [36, 37, 38].

This concept of adversarial attacks is potentially applicable to promoting
cooperation in the prisoner’s dilemma game. Even when defection is the
rational choice for individual players, strategically generated small
perturbations applied to the network may facilitate the evolution of
cooperation. This approach could offer a low-cost method to promote
cooperation without requiring large-scale network alterations.

Based on this idea, this study proposes a method to promote cooperation
evolution in the prisoner’s dilemma game in complex networks using adversarial
attacks. The method adds small, strategically generated perturbations to link
weights to foster cooperation, while also demonstrating the approach’s ability
to inhibit it. The proposed method is numerically evaluated across various
networks, including representative models and real-world social networks.
Additionally, the performance of this adversarial attack approach is compared
with an existing method for link weight adjustment, demonstrating its
superiority in promoting cooperation.

##  2 Methods

###  2.1 Prisoner’s dilemma game in complex networks

This study employs a prisoner’s dilemma game in a network of $N$ nodes
(players), based on settings used in previous research [1, 12, 13]. At time
$t$, each player adopts either a cooperative strategy (C;
$\boldsymbol{\sigma}_{i}(t)=(1~{}0)$) or a defective strategy (D;
$\boldsymbol{\sigma}_{i}(t)=(0~{}1)$), and gains payoffs by playing the game
with neighboring players. Following previous studies [19, 21, 23, 25, 26, 39,
40, 27] on the prisoner’s dilemma game in weighted networks, which focus on
adjusting link weights, the total payoff for player $i$ at time $t$ is
calculated as follows:

$P_{i}(t)=\sum_{j=1}^{N}A_{ij}\left[\boldsymbol{\sigma}_{i}(t)\boldsymbol{M}\boldsymbol{\sigma}_{j}^{\top}(t)\right]$
(1)

Here, $A_{ij}$ represents an element of the weighted adjacency matrix
$\boldsymbol{A}$ of the network, indicating player $j$’s influence on player
$i$. For simplicity, this study considers complex networks with bidirectional
links and uniform link weights of 1 (i.e., each player exerts equal influence
on each neighboring player). Thus, $A_{ij}=A_{ji}=1$ if a relationship exists
between players $i$ and $j$, and $A_{ij}=A_{ji}=0$ otherwise. Self-loops are
excluded (i.e., $A_{ii}=0$ for $i=1,\dots,N$). $\boldsymbol{M}$ is the game’s
payoff matrix. This study considers the weak prisoner’s dilemma game, which is
commonly used in previous research [4, 6, 12, 13, 27, 17, 41, 42, 43, 5, 44,
39, 20, 19, 15, 14, 45]. The payoff matrix is simplified as follows:

$\boldsymbol{M}=\left({\begin{array}[]{cc}1&0\\\ b&0\end{array}}\right),$

where $b$ represents the advantage of defectors over cooperators, constrained
by $1<b<2$ in the prisoner’s dilemma game [4, 6, 12, 13].

Player $i$ determines their strategy for the next time step ($t+1$) based on
their total payoff and those of their neighbors (strategy update).
Specifically, player $i$ imitates the strategy of player $j$, selected with a
probability proportional to the link weight ($A_{ij}/\sum_{h}A_{ih}$),
according to the following Fermi rule [46, 23, 26, 44]:

$\phi_{ij}=\frac{1}{1+\exp\left[\left\\{P_{i}(t)-P_{j}(t)\right\\}/K\right]}.$
(2)

Here, $K$ represents environmental noise in the strategy adoption process,
including irrationality and errors. Following previous studies, $K$ is set to
0.1 [5, 19, 26, 44, 20, 14]. It has been confirmed that using an alternative
update rule, such as the proportional rule employed in other studies [12, 13,
14, 25], does not alter the conclusions of this study.

Players’ strategies are updated synchronously, and this update is repeated
$t_{\max}$ times. It has been verified that asynchronous strategy updates do
not change the study’s conclusions. Previous research has shown that whether
strategy updates are synchronous or asynchronous does not qualitatively affect
the results [13, 47].

###  2.2 Adversarial attacks

Following the approach of adversarial attacks on voter model dynamics [36,
37], this study considers adversarial attacks that drives players’ strategies
towards a target state in the next time step (i.e., $t+1$). For computational
convenience, player $i$’s strategy at time $t$ is defined as
$x_{i}(t)=(+1~{}-1)\boldsymbol{\sigma}^{\top}_{i}(t)$. Thus,
$x_{i}(t)=+1~{}(-1)$ if player $i$ chooses strategy C (D) at time $t$. Let
$x_{i}^{*}=\\{+1,-1\\}$ be the target state for players, and the attack aims
to achieve $x_{i}(t+1)=x_{i}^{*}$ for $i=1,\dots,N$. The attacks are applied
by minimizing the energy $E=-N^{-1}\sum_{i=1}^{N}x^{*}_{i}x_{i}(t+1)$ (i.e.,
the negative value of the correlation coefficient between the observed state
and target state). To promote cooperation, $x^{*}_{i}=+1$ is set for
$i=1,\dots,N$. Consequently, the energy $E$ to be minimized becomes:

$E=-\frac{1}{N}\sum_{i=1}^{N}x_{i}(t+1).$

To minimize $E$, the study considers temporally varying link weights.
Specifically, the gradient descent method is used to add perturbations to
$\boldsymbol{A}$ at each time step. Assuming that link weights between
unconnected node pairs in the original network cannot be changed,
perturbations are added to the link weights of node pairs $i$ and $j$ where
$A_{ij}\neq 0$ at time $t$ as follows:

$A_{ij}^{*}(t)=A_{ij}-\epsilon\frac{\partial E}{\partial A_{ij}}.$

Here, $\epsilon$ is a small value. Note that setting $\epsilon>0$ results in
adversarial attacks to promote cooperation (targeting C), while $\epsilon<0$
leads to attacks inhibiting cooperation or, equivalently, promoting defection
(targeting D).

However, the gradient $\partial E/\partial A_{ij}$ cannot be directly
(analytically) obtained. Drawing from previous studies [36, 38], the following
mean-field time evolution of player $i$’s strategy in the prisoner’s dilemma
game is considered:

$x_{i}(t+1)=\sum_{j=1}^{N}\frac{A_{ij}}{\sum_{h=1}^{N}A_{ih}}\left[\phi_{ij}x_{j}(t)+(1-\phi_{ij})x_{i}(t)\right]$

In this case, the gradient becomes:

$\begin{split}\frac{\partial E}{\partial
A_{ij}}=\frac{1}{N\left(\sum_{h=1}^{N}A_{ih}\right)^{2}}&\left[A_{ij}\phi_{ij}\\{x_{i}(t)-x_{j}(t)\\}\sum_{\begin{subarray}{c}h=1\\\
h\neq j\end{subarray}}^{N}A_{ih}\right.\\\
&\left.\qquad-\sum_{\begin{subarray}{c}h=1\\\ h\neq
j\end{subarray}}^{N}A_{ih}\phi_{ih}\\{x_{i}(t)-x_{h}(t)\\}\right].\end{split}$

However, directly using this gradient for adversarial attacks is not practical
as the exact values of $\phi_{ij}$ are unknown. Therefore, following previous
research [36, 37], the study utilizes the optimal maximum norm constrained
perturbation (perturbation based on the sign of the gradient) [32]. Since
$A_{ij}\geq 0$, $x_{i}(t)\in\\{+1,-1\\}$, and $N>0$, the sign of the gradient
is preserved even if $\phi_{ij}$ is replaced with another $\phi^{\prime}_{ij}$
that monotonically increases with respect to $P_{i}(t)-P_{j}(t)$. Thus, using
such a $\phi^{\prime}_{ij}$, the sign of the gradient can be correctly
estimated. Specifically, perturbations are added to each element of
$\boldsymbol{A}$ based on the sign of the gradient as follows:

$A_{ij}^{\mathrm{adv}}(t)=A_{ij}-\epsilon\times\mathrm{sign}\left(\frac{\partial
E}{\partial A_{ij}}\right).$

Using the sign of the gradient in this way also has the advantage of
controlling the strength of the perturbations with $\epsilon$. Note that
$\boldsymbol{A}^{\mathrm{adv}}$ may be asymmetric.

The adversarial attacks on the prisoner’s dilemma game are primarily executed
by replacing $\boldsymbol{A}$ with $\boldsymbol{A}^{\mathrm{adv}}$ in Equation
(1). This change indirectly affects Equation (2) through the modified payoffs.

###  2.3 Comparison with other methods

To assess the efficacy of the proposed adversarial attacks, comparisons are
made with random attacks and Li et al.’s link weight adjustment method [25].

Random attacks serve as a control condition, providing a basis for evaluating
the significance of gradient consideration in the adversarial attacks. These
attacks involve adding perturbations to $\boldsymbol{A}$ according to the
following formula: $A_{ij}^{\mathrm{rnd}}(t)=A_{ij}+\epsilon\times s$ for all
pairs $(i,j)$ where $A_{ij}\neq 0$. In this context, $s$ represents a random
variable sampled from the set $\\{-1,+1\\}$. The implementation of the random
attacks primarily involves substituting $\boldsymbol{A}$ with
$\boldsymbol{A}^{\mathrm{rnd}}$ in Equation (1).

Li et al.’s method is employed as a representative link weight adjustment
technique for promoting cooperation. This selection is justified by its status
as a state-of-the-art approach and, analogous to the proposed method, its
capacity to regulate the magnitude of link weight modifications (perturbation
strength), thus facilitating a fair comparison. In Li et al.’s method, an
adaptive adjustment is implemented by comparing a player’s total payoff with
the mean total payoff of their neighboring players. The link weight is
increased (decreased) by a unit amount $d$ if the player’s total payoff
exceeds (falls below) the average. Specifically, this method initializes
$A_{ij}^{\mathrm{Li}}(0)\leftarrow A_{ij}$ and, assuming $|V(i)|>0$ for
$i=1,\dots,N$, where $V(i)$ denotes the set of the neighbors of node $i$,
iteratively updates link weights of node pairs $i$ and $j$ where $A_{ij}\neq
0$ according to the following equation:

$A_{ij}^{\mathrm{Li}}(t+1)\leftarrow\mathrm{clip}_{A{ij}\pm\epsilon}\left[A_{ij}^{\mathrm{Li}}(t)+d\times\mathrm{sign}\left(P_{i}(t)-\frac{\sum_{j\in
V(i)}P_{j}(t)}{|V(i)|}\right)\right].$ (3)

Here, $\mathrm{clip}_{A_{ij}\pm\epsilon}(x)$ is a function that constrains $x$
to the interval $[A_{ij}-\epsilon,A_{ij}+\epsilon]$. Note that
$0\leq\epsilon<A_{ij}$ (given the configuration of $\boldsymbol{A}$ in this
study, $0\leq\epsilon<1$). In accordance with previous research, $d$ is set to
0.01. The prisoner’s dilemma game based on this method is primarily executed
by substituting $\boldsymbol{A}$ with $\boldsymbol{A}^{\mathrm{Li}}$ in
Equation (1).

##  3 Results

To evaluate the effect of adversarial attacks on the prisoner’s dilemma game,
simulations were first conducted on complex networks generated from three
representative network models: the Erdős–Rényi (ER) model [7, 8],
Barabási–Albert (BA) [7, 10], and Watts-Strogatz (WS) models [9].

The ER model, widely used in network science, generates random networks by
establishing links between $L$ randomly selected node pairs from all possible
pairs. In this model, the node degree follows a Poisson distribution with mean
$\langle k\rangle=2L/N$. However, real-world social networks exhibit nonrandom
structures, characterized by highly clustered subnetworks and heterogeneous
(power-law-like) degree distributions [7, 10, 11]. Thus, BA and WS models were
considered. The BA model generates scale-free random networks by iteratively
adding new nodes and connecting them to $m$ existing nodes using a
preferential attachment mechanism. This process results in a degree
distribution $P(k)$ that follows a power law [$P(k)\propto k^{-3}$]. It should
be noted that $\langle k\rangle=2m$ for $N\gg 0$. The WS model generates
small-world networks by randomly rewiring links in a one-dimensional lattice,
where each node initially has $k$ ($=\langle k\rangle$) neighbors, with
probability $p_{\mathrm{WS}}$. The resulting networks exhibit higher
clustering coefficients than those expected from ER networks. In this study,
$p_{\mathrm{WS}}=0.05$, following [36].

In the simulations, unless otherwise specified, networks with $N=1000$ nodes
and an average degree of $\langle k\rangle=8$ were used. It has been confirmed
that using different values of $N$ or $\langle k\rangle$ does not
qualitatively affect the results. The initial state ($t=0$) consisted of an
equal number of cooperators and defectors, and the simulation time $t_{\max}$
was set to 1000 steps. The proportion of cooperators $\rho$ in the equilibrium
state was calculated by averaging over the last 100 steps. It has been
verified that this $t_{\max}$ setting is sufficient for $\rho$ to reach
equilibrium, and that moderate changes in the number of steps used for
averaging do not affect the results. The structure of model networks and the
initial distribution of cooperators and defectors within the networks depend
on random number generation. To minimize the impact of individual network
structures and initial states’ randomness, and to more accurately evaluate the
effect of adversarial attacks, the value of $\rho$ is reported as the average
of results obtained from 100 independent trials.

Figure 1 demonstrates that adversarial attacks can both promote and inhibit
cooperation across different network models (ER, BA, and WS). For promoting
cooperation ($\epsilon>0$), $\rho$ maintains higher values even for large $b$
compared to the case without attacks ($\epsilon=0$). This effect is observed
in all network types, with ER networks (Figure 1a) showing a substantial
increase in $\rho$ from $\sim 0.2$ ($\epsilon=0$) to $\sim 1.0$
($\epsilon=0.15$) at $b=1.25$. The range of b where cooperation dominates
($\rho>0.5$) is significantly broadened, as seen in the extension of the upper
limit of $b$ in ER networks from $\sim 1.15$ ($\epsilon=0$) to $\sim 1.5$
($\epsilon=0.2$).

Similar promotion of cooperation is observed in BA and WS networks (Figures 1b
and 1c), albeit with some differences due to their unique structural
properties. In BA networks, which inherently promote cooperation, the effect
is particularly pronounced. For example, at $b=1.95$, $\rho$ increases from
$\sim 0.2$ ($\epsilon=0$) to $\sim 0.9$ ($\epsilon=0.2$), with the upper limit
of $b$ for dominant cooperation extending beyond $1.95$ when $\epsilon\geq
0.1$. WS networks, despite their tendency to inhibit cooperation, also show
significant improvements. At $b=1.15$, $\rho$ increases from $\sim 0.05$
($\epsilon=0$) to $\sim 1.0$ ($\epsilon=0.15$), with the upper limit of $b$
for dominant cooperation increasing from $\sim 1.05$ ($\epsilon=0.0$) to $\sim
1.3$ ($\epsilon=0.2$).

Figure 1:  Line plots of the proportion of cooperators $\rho$ versus the
advantage of defectors over cooperators $b$ for different $\epsilon$ in (a)
ER, (b) BA, and (c) WS networks.

The effect of adversarial attacks is also remarkable for the inhibition of
cooperation ($\epsilon<0$) (Figure 1). While $\rho$ decreases with increasing
perturbation strength $|\epsilon|$ in adversarial attacks, the decrease due to
random attacks is limited across all network types. In ER networks, at
$b=1.1$, $\rho$ decreases from $\sim 1.0$ ($\epsilon=0$) to $\sim 0.1$
($\epsilon=-0.1$), with the upper limit of $b$ for dominant cooperation
falling below $1.05$ when $\epsilon=-0.1$. BA networks, despite their inherent
cooperation-promoting nature, show a significant reduction in cooperation
under adversarial attacks. For instance, at $b=1.5$, $\rho$ decreases from
$\sim 0.9$ ($\epsilon=0$) to $\sim 0.05$ ($\epsilon=-0.2$), with the upper
limit of $b$ for dominant cooperation decreasing to $\sim 1.2$. In WS
networks, the effect is even more pronounced, with $\rho$ dropping to $\sim
0.0$ at $b=1.05$ when $\epsilon\leq-0.1$.

The superiority of adversarial attacks becomes evident when compared to random
attacks and Li et al.’s method (Figure 2). For promoting cooperation
($\epsilon>0$), adversarial attacks show a rapid increase in $\rho$ with
increasing $\epsilon$, while random attacks fail to increase $\rho$ in any of
the model networks examined. Li et al.’s method shows limited effectiveness,
with slight increases in $\rho$ for ER and BA networks (up to $\sim 0.15$ in
ER networks) but no increase in WS networks.

Figure 2:  Line plots of the proportion of cooperators $\rho$ versus
perturbation strength $\epsilon$ for adversarial attacks promoting cooperation
($\epsilon>0$) in (a) ER, (b) BA, and (c) WS networks. The advantage of
defectors over cooperators $b$ is 1.35, 1.95, and 1.2 for (a), (b), and (c),
respectively.

The trend of $\rho$ increase due to adversarial attacks varies among network
structures, with BA networks showing an almost linear increase, while ER and
WS networks exhibit a delayed increase, particularly pronounced in WS
networks. This is likely due to the cooperation-suppressing effect of WS
networks compared to other networks, as well as the $\epsilon$ being too
small. Indeed, in previous studies, relatively large $\epsilon$ values were
set to sufficiently promote cooperation. Note that in Li et al.’s method, link
weights are incrementally increased or decreased, but it has been confirmed
that the change in link weights reaches $\epsilon$ in early time steps.

These results suggest that adversarial attacks can promote cooperation more
efficiently. However, note that the trend of $\rho$ increase due to
adversarial attacks varies somewhat depending on the network structure. For
example, in BA networks, the increase in $\rho$ is almost linear with respect
to the increase in $\epsilon$, whereas in ER and WS networks, there is a delay
in this increase. This delay trend is particularly pronounced in WS networks.
These differences are likely attributable to the inherent characteristics of
each network model (as mentioned earlier, BA networks have an intrinsic effect
of promoting cooperation, while WS networks have an effect of suppressing
cooperation).

For inhibiting cooperation ($\epsilon<0$) (Figure 3), adversarial attacks
consistently and effectively suppress cooperation across all model networks.
In contrast, random attacks show limited or even counterproductive effects,
particularly in WS networks where they tend to increase $\rho$, likely due to
the weakening of cooperation-inhibiting shortcuts. The inhibitory effect of
adversarial attacks, while consistent, is slightly influenced by network
structure. ER and WS networks show rapid decreases in $\rho$ with increasing
$|\epsilon|$, whereas the decrease is more gradual in BA networks, again
reflecting their intrinsic cooperation-promoting nature.

Figure 3:  Line plots of the proportion of cooperators $\rho$ versus
perturbation strength $|\epsilon|$ for adversarial attacks inhibiting
cooperation ($\epsilon<0$) in (a) ER, (b) BA, and (c) WS networks. The
advantage of defectors over cooperators $b$ is 1.1, 1.5, and 1.05 for (a),
(b), and (c), respectively.

The effect of adversarial attacks is further evaluated using real-world social
networks (Figure 4): Facebook [48], Advogato [49, 50], AnyBeat [49, 51], and
HAMSTERster networks [49]. These networks reflect actual human relationships
obtained from online social media platforms and web services, and possess
diverse structural characteristics (e.g., high degree heterogeneity and high
clustering coefficients; see Table 1 in [37] for details). Note that these
networks are undirected. For simplicity, simulations used the largest
connected component in each network, with all link weights set to 1. The
effect of adversarial attacks on real-world networks was found to be
qualitatively consistent with the results from the aforementioned model
networks. The results are qualitatively consistent with those from model
networks, demonstrating both promotion ($\epsilon>0$) and inhibition
($\epsilon<0$) of cooperation. However, the AnyBeat network shows
comparatively smaller effects, highlighting a limitation of the proposed
method. This is likely due to its sparse nature and the method’s constraint of
modifying only existing link weights. In this adversarial attack approach,
perturbations are limited to node pairs already connected in the original
network, potentially reducing the method’s effectiveness in sparse networks.

In the adversarial attack method of this study, only the link weights of node
pairs connected in the original network are modified, thus limiting the effect
of attacks in sparse networks.

Figure 4:  Line plots of the proportion of cooperators $\rho$ versus the
advantage of defectors over cooperators $b$ for different $\epsilon$ in (a)
Facebook ($N=4039$ and $\langle k\rangle=43.7$), (b) Advogato ($N=5054$ and
$\langle k\rangle=16.6$), (c) AnyBeat ($N=12645$ and $\langle k\rangle=7.8$),
and (d) HAMSTERster networks ($N=2000$ and $\langle k\rangle=16.1$).

##  4 Discussion

This study investigated the impact of strategic manipulation of link weights,
applying the concept of adversarial attacks, on the evolution of cooperation
in the prisoner’s dilemma game in complex networks. It revealed that even
slight perturbations can effectively promote or inhibit cooperation (Figures
1–4). This effect was confirmed not only in model networks (Figure 1) but also
in real-world social networks (Figure 4), demonstrating the versatility and
robustness of adversarial attacks. Conventional approaches, especially link
weight adjustment methods represented by Li et al.’s technique, required
relatively large changes in link weights to promote cooperation. In contrast,
the approach utilizing the concept of adversarial attacks proposed in this
study enables effective evolution of cooperation with minimal perturbations
(Figure 2). This point is particularly crucial because, in real-world social
networks, the structure and link weights are determined by individuals’
intrinsic properties and relationships, making large changes unrealistic and
only small modifications permissible.

In connection with this, this study suggests that large heterogeneity in link
weights [19, 20, 21, 22, 23, 27] is not necessary to promote cooperation. The
proposed adversarial attack method demonstrates that cooperation can be
effectively promoted with minimal weight modifications. This provides a new
perspective on link weights’ role in cooperation evolution, suggesting that
the ability to fine-tune network interactions may be more crucial than static
weight heterogeneity.

Link weights correspond to factors such as frequency of contact or level of
trust between individuals [21, 24, 19, 23]. Making slight modifications to
these is relatively simple in social networking services (SNS), for instance.
Operators with access to social network data and user behavior patterns can
potentially promote cooperative behavior across the network by manipulating
the frequency of other users’ posts appearing in individual timelines or
altering scores related to personal trust levels. Notably, these link weight
modifications can be relatively easily implemented as they are based solely on
individual behavior patterns and payoff information (Equation (2.2)).

However, caution is necessary regarding the effectiveness of this adversarial
attack in inhibiting cooperation (promoting defection; Figure 3). Operators or
malicious third parties with access to social network data and user behavior
patterns could potentially suppress the evolution of cooperation and cause
social disruption without users’ knowledge. This point suggests a new
vulnerability in digital society. Furthermore, the application of this method
extends beyond social media to various social interaction settings, such as
organizational communication tools and online learning platforms. This
presents both opportunities for promoting cooperative behavior and ethical
challenges concerning individual autonomy. These findings have important
implications for the design and operation of social systems, highlighting the
need for careful ethical consideration and appropriate regulation in the use
of these technologies.

The adversarial attack concept is expected to continue playing a crucial role
in developing approaches to promote cooperation. Previous methods for
promoting the evolution of cooperation have primarily utilized either payoffs
(e.g., [25, 26, 24]) or individual strategies (cooperation or defection; e.g.,
[27, 29, 30, 28]), and have often been designed empirically. By introducing
the adversarial attack concept, it becomes possible to reframe this problem as
a kind of optimization problem, enabling a more systematic approach
derivation. Indeed, the method demonstrated in this study achieves more
efficient promotion of cooperation by combining both payoffs and individual
strategies. This suggests that the framework of adversarial attacks possesses
the flexibility to comprehensively consider multiple factors.

Further development of this approach may lead to more efficient and
sophisticated methods for promoting cooperation. For instance, the current
study has constraints such as approximating the prisoner’s dilemma game with
mean-field time evolution to estimate gradients, and only being able to
estimate the sign of the gradient. Due to these constraints, while the
proposed method was effective for representative model networks and several
real-world networks, its effect may be limited for specific network structures
or under dynamic conditions. Future research directions include developing
more accurate gradient estimation methods. The exploration of more precise
gradient estimation approaches that go beyond mean-field approximation is
anticipated. Moreover, developing sparse attack strategies that work
effectively with manipulation of fewer links is another important challenge.
This could lead to reduced computational costs and easier implementation.

Furthermore, the adversarial attacks examined in this study are limited to
complex networks where relationships between individuals are bidirectional and
all link weights are equal. Therefore, investigating adversarial attacks in
complex networks where interpersonal relationships are asymmetric [18, 52] and
link weights vary significantly [19, 20, 21, 22, 23] is of great interest from
the perspective of application to more realistic social network models. In
this regard, devising adaptive attack methods that can respond to temporal
networks [29, 45] and multiplex networks [53, 54, 55] is also an important
research challenge, as real-world social networks are not static but
constantly changing.

Applying these methods to different game-theoretic models [1, 2] (e.g., public
goods games [56, 57] or coordination games [58, 59]) or situations where
multiple games coexist [60, 61]. This could deepen our understanding of
cooperation promotion and inhibition mechanisms in various social dilemma
situations. Additionally, it is important to verify these methods in
situations with players having diverse strategies [56, 62]. In real societies,
various behavioral strategies exist beyond simple cooperation and defection,
so examining the effects of adversarial attacks under such complex conditions
is expected to yield more realistic insights.

Moreover, verification in social networks composed of actual humans is
essential [16, 63]. By examining the extent to which findings from simulations
and theoretical models are applicable to real human behavior, the practicality
and social impact of this research can be evaluated.

Advancing research on adversarial attacks in evolutionary games has the
potential to enhance our understanding of cooperation and competition dynamics
in complex social systems, contributing to more effective social system
design.

##  5 Conclusions

This study introduces a novel approach to modulating cooperation in complex
networks through adversarial attacks on the prisoner’s dilemma game. We
demonstrate that minimal, strategic perturbations can significantly influence
cooperation evolution across diverse network structures, challenging the
necessity of large link weight heterogeneity. The method’s efficacy in both
theoretical and real-world networks highlights its potential for practical
application. By framing cooperation promotion as an optimization problem, we
offer a more systematic approach with implications for social system design.
Future research should explore efficient strategies for varied network
topologies, extend to other evolutionary games, and validate results
empirically. This work advances our understanding of cooperation dynamics and
provides a tool for shaping collective behavior, while emphasizing the
importance of ethical considerations in digital social systems.

## Data availability statement

The data and code that support the results presented in this study are openly
available in the author’s GitHub repository (github.com/kztakemoto/AdvGame).

This study was supported by JSPS KAKENHI (Grant Number 21H03545).

## References

## References

  * [1] Szabó G and Fath G 2007 Physics reports 446 97–216
  * [2] Perc M, Jordan J J, Rand D G, Wang Z, Boccaletti S and Szolnoki A 2017 Physics Reports 687 1–51
  * [3] Ohtsuki H and Nowak M A 2006 Journal of theoretical biology 243 86–97
  * [4] Nowak M A and May R M 1992 Nature 359 826–829
  * [5] Szabó G and Tőke C 1998 Physical Review E 58 69
  * [6] Nowak M A and May R M 1993 International Journal of bifurcation and chaos 3 35–78
  * [7] Albert R and Barabási A L 2002 Reviews of modern physics 74 47
  * [8] Takemoto K and Oosawa C 2012 Introduction to Complex Networks: Measures, Statistical Properties, and Models (John Wiley & Sons, Ltd) chap 2, pp 45–75 ISBN 9781118346990
  * [9] Watts D J and Strogatz S H 1998 Nature 393 440–442
  * [10] Barabási A L and Albert R 1999 Science 286 509–512
  * [11] Barabási A L and Bonabeau E 2003 Scientific American 288 60–69
  * [12] Santos F C and Pacheco J M 2005 Physical review letters 95 098104
  * [13] Santos F C, Rodrigues J and Pacheco J M 2006 Proceedings of the Royal Society B: Biological Sciences 273 51–55
  * [14] Wu Z X, Guan J Y, Xu X J and Wang Y H 2007 Physica A: Statistical Mechanics and its Applications 379 672–680
  * [15] Szolnoki A, Perc M and Danku Z 2008 Physica A: Statistical Mechanics and its Applications 387 2075–2082
  * [16] Gracia-Lázaro C, Ferrer A, Ruiz G, Tarancón A, Cuesta J A, Sánchez A and Moreno Y 2012 Proceedings of the National Academy of Sciences 109 12922–12926
  * [17] Abramson G and Kuperman M 2001 Physical Review E 63 030901
  * [18] Kim B J, Trusina A, Holme P, Minnhagen P, Chung J S and Choi M 2002 Physical Review E 66 021907
  * [19] Du W B, Zheng H R and Hu M B 2008 Physica A: Statistical Mechanics and its Applications 387 3796–3800
  * [20] Meng X, Sun S, Li X, Wang L, Xia C and Sun J 2016 Physica A: Statistical Mechanics and its Applications 442 388–396
  * [21] Iwata M and Akiyama E 2016 Physica A: Statistical Mechanics and its Applications 448 224–234
  * [22] Hong-Yan C, Hai-Hong L, Qiong-Lin D and Jun-Zhong Y 2010 Chinese Physics Letters 27 070201
  * [23] Buesser P and Tomassini M 2012 Physical Review E 85 016107
  * [24] Cao L, Ohtsuki H, Wang B and Aihara K 2011 Journal of theoretical biology 272 8–15
  * [25] Li X, Sun S and Xia C 2019 Applied Mathematics and Computation 361 810–820
  * [26] Liu C, Wang J, Li X and Xia C 2020 Physica A: Statistical Mechanics and its Applications 554 124691
  * [27] Li D, Ma J, Han D, Sun M, Tian L and Stanley H E 2017 Scientific reports 7 45237
  * [28] Santos F C, Pacheco J M and Lenaerts T 2006 PLoS computational biology 2 e140
  * [29] Fu F, Wu T and Wang L 2009 Physical Review E 79 036101
  * [30] Lee H W, Malik N and Mucha P J 2018 Journal of complex networks 6 1–23
  * [31] Szegedy C, Zaremba W, Sutskever I, Bruna J, Erhan D, Goodfellow I J and Fergus R 2014 Intriguing properties of neural networks 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings ed Bengio Y and LeCun Y URL http://arxiv.org/abs/1312.6199
  * [32] Goodfellow I J, Shlens J and Szegedy C 2015 Explaining and harnessing adversarial examples 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings ed Bengio Y and LeCun Y URL http://arxiv.org/abs/1412.6572
  * [33] Yuan X, He P, Zhu Q and Li X 2019 IEEE Transactions on Neural Networks and Learning Systems 30 2805–2824 ISSN 2162-237X URL https://ieeexplore.ieee.org/document/8611298/
  * [34] Hirano H, Minagi A and Takemoto K 2021 BMC medical imaging 21 9
  * [35] Hirano H, Koga K and Takemoto K 2020 Plos one 15 e0243963
  * [36] Chiyomaru K and Takemoto K 2022 Physical Review E 106 014301
  * [37] Chiyomaru K and Takemoto K 2023 Journal of Physics: Complexity 4 025009
  * [38] Mizutaka S 2023 Journal of Physics: Complexity 4 035009
  * [39] Wen-Bo D, Xian-Bin C, Lin Z and Hong Z 2009 Chinese Physics Letters 26 058701
  * [40] Zhai C, Zhang H T, Zhao Y, Chen M Z, Rong Z H and Wang B H 2010 Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 81 036112
  * [41] Wang X P, Jiang L L and Wang B H 2012 Chinese Physics B 21 070210
  * [42] Devlin S and Treloar T 2009 Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 80 026105
  * [43] Szabó G and Szolnoki A 2009 Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 79 016106
  * [44] Liu S, Zhang L and Wang B 2020 New Journal of Physics 22 013034
  * [45] Li A, Zhou L, Su Q, Cornelius S P, Liu Y Y, Wang L and Levin S A 2020 Nature communications 11 2259
  * [46] Zhi-Xi W, Xin-Jian X and Ying-Hai W 2006 Chinese Physics Letters 23 531
  * [47] Hauert C and Doebeli M 2004 Nature 428 643–646
  * [48] Leskovec J and Mcauley J 2012 Advances in neural information processing systems 25
  * [49] Rossi R and Ahmed N 2015 The network data repository with interactive graph analytics and visualization Proceedings of the AAAI conference on artificial intelligence vol 29
  * [50] Massa P, Salvetti M and Tomasoni D 2009 Bowling alone and trust decline in social network sites 2009 Eighth IEEE International Conference on Dependable, Autonomic and Secure Computing (IEEE) pp 658–663
  * [51] Fire M, Puzis R and Elovici Y 2013 Handbook of computational approaches to counterterrorism 283–300
  * [52] Gao L, Li Y, Wang Z and Wang R W 2022 Applied Mathematics and Computation 412 126590
  * [53] Gómez-Gardenes J, Reinares I, Arenas A and Floría L M 2012 Scientific reports 2 620
  * [54] Takesue H 2021 Applied Mathematics and Computation 388 125543
  * [55] Wang Z, Wang L, Szolnoki A and Perc M 2015 The European physical journal B 88 1–15
  * [56] Santos F C, Santos M D and Pacheco J M 2008 Nature 454 213–216
  * [57] Gómez-Gardenes J, Vilone D and Sánchez A 2011 Europhysics Letters 95 68003
  * [58] Broere J, Buskens V, Weesie J and Stoof H 2017 Scientific reports 7 17016
  * [59] Raducha T and San Miguel M 2023 Scientific Reports 13 11818
  * [60] Szolnoki A and Perc M 2014 Europhysics Letters 108 28004
  * [61] Venkateswaran V R and Gokhale C S 2019 Proceedings of the Royal Society B 286 20190900
  * [62] Nowak M and Sigmund K 1993 Nature 364 56–58
  * [63] Rand D G, Arbesman S and Christakis N A 2011 Proceedings of the National Academy of Sciences 108 19193–19198


# Combinatorial optimization solving by coherent Ising machines based on
spiking neural networks

Bo Lu  School of Artificial Intelligence, Beijing Normal University, Beijing
100875, China  Yong-Pan Gao  School of Electronics Engineering, Beijing
University of Posts and Telecommunications, Beijing 100876, China  Kai Wen
Beijing QBoson Quantum Technology Co., Ltd., Beijing 100015, China  Chuan Wang
<EMAIL_ADDRESS>School of Artificial Intelligence, Beijing Normal
University, Beijing 100875, China

###### Abstract

Spiking neural network is a kind of neuromorphic computing which is believed
to improve on the level of intelligence and provide advabtages for quantum
computing. In this work, we address this issue by designing an optical spiking
neural network and prove that it can be used to accelerate the speed of
computation, especially on the combinatorial optimization problems. Here the
spiking neural network is constructed by the antisymmetrically coupled
degenerate optical parametric oscillator pulses and dissipative pulses. A
nonlinear transfer function is chosen to mitigate amplitude inhomogeneities
and destabilize the resulting local minima according to the dynamical behavior
of spiking neurons. It is numerically proved that the spiking neural network-
coherent Ising machines has excellent performance on combinatorial
optimization problems, for which is expected to offer a new applications for
neural computing and optical computing.

##  I Introduction

Combinatorial optimization problems are of great importance in various fields
of computer science. During the past decades, they are ubiquitous in various
field, such as the drug design Kitchen  _et al._ (2004), financial management
Soler-Dominguez  _et al._ (2017), and circuit design Barahona  _et al._
(1988). Due to the rapid growth of the computing power, the computation time
of these problems on conventional computers tends to increase exponentially
with the size of the problem due to the explosion in the number of
combinations Arora and Barak (2009). The problem in the complexity class NP
could be formulated as an Ising problem with only polynomial overhead, and the
loss functions can be mapped to the Ising model to complete the search from
the initial state to the ground state, which has been implemented by several
physical systems such as quantum annealing Kadowaki and Nishimori (1998); Das
and Chakrabarti (2008), optical oscillators Wang  _et al._ (2013),
nanomagnetic-net arrays Sutton  _et al._ (2017); Saccone  _et al._ (2022).

Among them, oscillator-based coherent Ising machines (CIM) are widely studied
due to the advantages Hamerly  _et al._ (2019) of high speed, programmability,
and parallel search. Each Ising spin in the CIM can be encoded as the phase
$\phi_{i}$ of the light pulse in an optical mode. Employing the degenerate
optical parametric oscillation, the phase values are enforced as the binary
spin values as $\phi_{i}=0,\pi$. During computation, the optical gain yields
oscillation of the pluses either in-phase or out-of-phase with respective to
the pump light. Meanwhile, the oscillators are coupled with tunable coupling
constant to form a network. The oscillation network can be implemented by
various ways, such as the optical oscillators Marandi  _et al._ (2014);
Inagaki  _et al._ (2016a); McMahon  _et al._ (2016); Inagaki  _et al._
(2016b); Honjo  _et al._ (2021), opto-electronic oscillators Fabian  _et al._
(2019); Cen  _et al._ (2020), spatial phase modulators Pierangeli  _et al._
(2019); Fang  _et al._ (2021); Sun  _et al._ (2022), and integrated photonic
chips Okawachi  _et al._ (2020).

As a quantum heuristic algorithm, CIM has excellent performance on small-scale
combinatorial optimization problems Hamerly  _et al._ (2019); McMahon  _et
al._ (2016); Wang  _et al._ (2013), but on large-scale problems, an important
technicality would arise in the model which is that the amplitude of the
oscillators are not equal. This phenomenon is usually call the amplitude
heterogeneity which prevents CIM from being correctly mapped to the Ising
model Leleu  _et al._ (2017), as shown in Fig.1(a). The spin configuration of
the ground state of the Ising model differs from that of the CIM. During the
calculation of CIM, the energy will be frozen in the excited state of the
Ising model, as shown in Fig.1(b). To address the problem of amplitude
heterogeneity, there are several method to fix these problem by adding a
feedback mechanism to reduce the unequality of the amplitude, such as the
additional uncorrelated driving signals Leleu  _et al._ (2017), nonlinear
functions Böhm  _et al._ (2021); Reifenstein  _et al._ (2021) and correcting
the amplitude by the error signal Leleu  _et al._ (2019); Kako  _et al._
(2020); Leleu  _et al._ (2021).

In addition, we turn our attention to the spiking neural network (SNN) which
has been simulated on the CIM based on the DOPOs neural network, and the
collective behavior of neural clusters has been studied Inagaki  _et al._
(2021). The spiking neurons which are a special type of neurons which performs
the signal processing in the brain. They are considered as the next generation
of neural architectures, and have been extensively studied in artificial
intelligence recently Maass (1997); Brette  _et al._ (2007); Izhikevich
(2003). For example, it can imitate the biological nervous system to produce
more intelligent systems, such as brain-like learning Mnih  _et al._ (2015);
Yu  _et al._ (2012).

In this work, we propose a new architecture of CIM, based on a spiking neuron
network composed of DOPO pulses and dissipative pulses. We find that it is
able to achieve the amplitude heterogeneity and also freeze in local minima
during computating for combinatorial optimization problems. By controlling the
dissipative parameters of spiking neurons composed of antisymmetrically
coupled DOPO and dissipative pulses, the rate and direction of contraction of
the phase space volumes can be controlled to bring the system out of the
frozen state. More generally, we prove the existence of neurons are similar to
the dynamic behavior of class-II neurons. A SNN-CIM algorithm is proposed to
achieve the instability of the local minimum by adjusting the parameters, and
the performance is compared with the algorithm that only uses the nonlinear
filter function. The results show that SNN-CIM has excellent performance and
has more potential on complex graph structures.

##  II Results

###  II.1 SNN-CIM model

Figure 1: Schematic diagram of SNN-CIM. a. Variation of energies with spin
configuration for homogeneous and inhomogeneous spin amplitudes. The red line
is the pump intensity, which does not excite the correct ground state spin
configuration. b. The time evolution of the energy, due to the amplitude
heterogeneity, is frozen in the excited state. The dashed line is the ground
state energy. c. Optical implementation of SNN-CIM. The pump is turned on to
generate type-x DOPO, and when it is turned off, the modulation generates
type-k pulses. d. Spiking neural network. IM: intensity modulator, PM: phase
modulator, SHG: second harmonic generation, PSA: phase sensitive
amplification, FPGA: field programmable gate array.

The model of SNN-CIM we use can be summarized in Fig.1(c). The system with
correction of amplitude heterogeneity can be implemented using a CIM with a
measurement-and-feedback (MFB) structure composed of DOPO pulses. The qubits
for computing are represented by the optical SNN pulses which propagate
cyclically in the fiber-ring cavity. For each SNN pulses, it can be divided
into two types of pulses by controlling the pump, one is generated by the
optical parametric oscillation with nonlinear gain process. It corresponds to
the case where the pump is turned on, denoted as $x$, and the other pulse is
the injected optical pulse with a dissipative process only, corresponding to
the case where the pump is turned off, denoted as $k$.

A spiking neuron consists of both a type-$x$ pulse and a type-$k$ pulse
coupled to each other with coupling strengths of $J_{xk}$ and $J_{kx}$
$(J_{kx}=-J_{xk})$, respectively. The spiking neurons can be connectted to
each other, as shown in Fig.1(d). Based on this principle, the dynamical
behavior of the $i$th neuron can be expressed as:

$\begin{split}&\frac{dx_{i}}{dt}=\alpha x_{i}-x_{i}^{3}+J_{xk}\beta
k_{i}+I_{ext},\\\ &\frac{dk_{i}}{dt}=-\beta k_{i}+J_{kx}x_{i},\\\
&I_{ext}=tanh(\epsilon\Sigma_{j}J_{ij}x_{j}).\end{split}$ (1)

where $x_{i}$ is in-phase components of DOPO amplitudes, $(k_{i})$ is
dissipative pulse amplitude. $\alpha$ denotes the nonlinear gain generated by
the pump,and $\beta$ is the inherent dissipation of the system. $J_{ij}$
represents the connection between the $i$th and the $j$th neuron, and
$\epsilon$ is the strenght of interaction.

For the case of the power above the threshold, the DOPO pulses undergo
spontaneous symmetry breaking. Then, the potential energies of the spiking
neurons can be represented by the positive and negative amplitudes of the
type-$x$ pulses. According to the constraints of the combinatorial
optimization problems, the connection of the neurons in the neural network can
be regarded as an Ising model that evolves a spin configuration into the
ground state.

###  II.2  The dynamics of the spiking neurons

Figure 2: Dynamic behavior of spiking neurons. a. The effect of pump strength
$\alpha$ and $I_{ext}$ on the stability of the fixed point. $\beta=0.3$. b.
and c. The time evolution and phase space trajectories of the type-$x$ DOPO
pulse and type-$k$ pulse at point m in (a), $\alpha=1$ and $I_{ext}=0.05$. d.
and e. Corresponding to n points, $\alpha=1$ and $I_{ext}=0.25$. f. The time
evolution of the type-$x$ DOPO pulse and type-$k$ pulse with the increase of
$I_{ext}$.

According to Eq.(1), we find the dynamic properties of a spiking neuron
depends on the linear stability of its fixed point. And the fixed points under
steady-state solution of Eq.(1) should satisfy the relation as:

$\begin{split}&k_{i}=\frac{1}{\beta}x_{i},\\\
&0=(\alpha-1)x_{i}-x_{i}^{3}+I_{ext}.\end{split}$ (2)

The solution satisfying Eq.(2) is denoted as $x_{0}$ and the linear stability
of a fixed point $(x_{0},x_{0}/\beta)$ can be obtained by the eigenvalues of
its Jacobian matrix:

$J_{i}=\left[\begin{array}[]{cc}\alpha-3x_{0}^{2}&-\beta\\\ 1&-\beta\\\
\end{array}\right].$ (3)

Here the eigenvalues $\lambda_{i}^{\pm}$ of the Jacobian matrix are given as
follows:

$\lambda_{i}^{\pm}=\frac{1}{2}(T\pm\sqrt{T^{2}-4\Delta}),$ (4)

where $\Delta$ and $T$ are the determinant and trace of the Jacobian matrix,
respectively. Then, we have
$\Delta=J_{xx}J_{kk}-J_{xk}J_{kx}=(3x_{0}^{2}-\alpha+1)\beta$ and
$T=J_{xx}+J_{kk}=\alpha-3x_{0}^{2}-\beta$. Numerical simulations give the
fixed-point linear $((x_{0},x_{0}/\beta))$ stability as a function of
$I_{ext}$ and the pump strength $\alpha$ when $\beta=0.3$, as shown in
Fig.2(a). The stability of those fixed points depends on
$Re(\lambda_{i}^{+})$, these fixed points become unstable when
$Re(\lambda_{i}^{+})>0$. However, the eigenvalues would become complex
conjugated when $T^{2}-4\Delta<0$. Hence, we mark $Re(\lambda_{i}^{+})=0$ and
$T^{2}-4\Delta=0$ in Fig.2(a), respectively. Then, the time evolution and the
oscillation phase space of those fixed points in the regions 2 and 3 are shown
in Fig.2(b)(c) and Fig.2(d)(e), respectively. Under a small disturbance of
$I_{ext}$, the system reciprocates on a limit cycle. And the system will
shrink to a stable focus under a large $I_{ext}$.

The spiking behavior of single neurons can be achieved by controlling the
parameter $I_{ext}$. When $I_{ext}$ increases, the amplitude of neuron pulses
gradually decreases. However, for the increment of $I_{ext}$ to a certain
value, the neuron pulses will no longer appear, as shown in Fig.2(f). This
process is consistant with the evolution direction of the arrow indicated by
Fig.2(a), that is, the system appears Andronov-Hopf (AH) bifurcation with the
increase of $I_{ext}$, which make the fixed points change from an unstable
limit cycle to a stable focused point. This proves that the optical pulse
network has the properties of the class-II neurons by antisymmetrically
coupling a pair of nonlinear gain DOPO pluse and dissipation optical pluse.

In addition, as shown in Fig.2(g), we find a smaller $I_{ext}$ could exhibit a
larger firing rate by taking a Fourier transform of the temporal evolution of
neurons under different $I_{ext}$. For the Ising model with antiferromagnetic
coupling, the evolution of the energy
$E_{Ising}=-\Sigma_{ij}\sigma_{i}\sigma_{j}/2$ is accompanied by an
amplification of the coupling term. In this model, $I_{ext}$ corresponds to
the coupling term on the $i$th spin, that is, the firing rate decreases
($I_{ext}$ increases) during the evolution of the Ising model. Moreover, the
cluster neurons can also be synchronized to a low firing rate.

###  II.3  Instability of the local Minimum and SNN-CIM Algorithms

Figure 3: SNN-CIM solves the MAX-CUT problem of G1. a. The effect of $\beta$
and $I_{ext}$ on the stability of the fixed point. $\alpha=1.0$. b. Evolution
of type-$x$ DOPO pulse amplitude with calculation time. c. Evolution of
type-$k$ pulse amplitude with calculation time. d. Energy and number of cuts
as a function of computation time.

Similar to the analysis of the linear stability of pump and $I_{ext}$ on the
fixed points, Fig.3(a) shows the effect of dissipation $\beta$ and $I_{ext}$
on the linear stability of those fixed points. Those fixed points that freeze
at local minima tend to be located in region 3. The stability of the local
minima can be influenced to escape from the frozen local minima by adjusting
the dissipation and coupling strength. To overcome the local minima and make
the system to escape from the frozen state, we choose $Re(\lambda_{i}^{+})>0$
by reducing the dissipation $\beta$ and coupling strength $\epsilon$ which
forces these fixed points into region 2 to destabilized. And then, the system
continues to search for the ground state with lower energy by gradually
increasing the coupling strength.

To give an example of the algorithm, we apply the SNN-CIM algorithm on the
MAX-CUT problem of G1.

In Fig.3(b), the time evolution of the type-$x$ DOPO and type-$k$ pulses, the
energy and the number of cuts of the system are given. Specifically, if the
energy is frozen in a local minima state, we tend to reduce $\beta$ which
makes the system unstable. Then, there is fluctuation of the energy and the
system continues to search for the lower energy state. Finally, if the system
energy is higher than the lowest energy before this moment, as in Fig. 3, the
optimal solutions could be obtained as $edges=11624$ near T=1350, after that
the system will be unstable until another lower energy is searched.

###  II.4  The performance of SNN-CIM Algorithm on MAX-CUT Problem

Here in this part, we describe the performance of SNN-CIM algorithm on the
Max-Cut problem, and then compare with other types of computing devices such
as those based on simulated annealing and conventional CIM. We discuss the
performance of these investigated devices on different size of the Ising
problems.

We choose the simulated annealing (SA) algorithm and the SNN-CIM algorithm on
the G1-G21 graph of 800 nodes, and the G22-G42 graph of 2000 nodes, and also
compared the performance with the conventional CIM and the CIM with signmoid
function optimization in Ref.Böhm  _et al._ (2021).

Figure 4: Computational performance of the SNN-CIM algorithm on the MAX-CUT
problem for a given problem. a. Relative error of the optimal value of
simulated annealing algorithm, CIM, CIM with sigmoid function, and SNN-CIM for
50 runs on the MAX-CUT problem of a graph with 800 nodes. b. More complex 2000
node situation.

The performance of the computational results are shown in Fig.4. Here we
choose the relative error $\Delta C=100(1-C/C_{opt})$ of the optimal value
during 50 runs as the result, and $C_{opt})$ denotes the best-known value
reported in literature. In the 800 nodes problem, CIM with sigmoid function
and SNN-CIM often find the solutions very close to the optimal value, while
the performance of the SA and CIM are relatively poorly. And, the mean values
of the relative errors of SNN-CIM and CIM with sigmoid function in the optimal
value of G1-G21 are $0.3\%$ and $0.4\%$, respectively. The performance of SNN-
CIM is slightly better than that of CIM with sigmoid function. Under the more
complex 2000-node problem, SNN-CIM is significantly closer to the optimal
value than the other three algorithms. In addition, for more complex problems
with random connectivity and non-uniform node degrees (G39,G40,G41,G42), SNN-
CIM with unstable local optima can also achieve better computational
performance. This shows that SNN-CIM will have more potential in more complex
combinatorial optimization problems.

Table 1:  Computational performance of SA, CIM, and SNN-CIM on GSET. TTS is time-to-solution of reach 98 percent of the best solution. The optimal solutions for SA, CIM and SNN-CIM are obtained under 10, 50 and 50 runs, respectively. Values in parentheses indicate that the results of $TTS_{98}$ cannot be obtained, and are replaced by $TTS_{95}$. The best maximum cuts known from Ref.Ma and Hao (2017) ID |  N |  $J_{ij}$ |  Edge |  $C_{opt}$ |  $C_{SA}$ |  $C_{CIM}$ |  $C_{SNN-CIM}$ |  $TTS_{SA}(s)$ |  $TTS_{CIM}(s)$ |  $TTS_{SNN-CIM}(s)$
---|---|---|---|---|---|---|---|---|---|---
1 |  800 |  $\\{0,1\\}$ |  19176 |  11624 |  11604 |  11597 |  11624 |  16.1514 |  3.8065 |  0.0373
6 |  800 |  $\\{0,\pm 1\\}$ |  19176 |  2178 |  2138 |  2149 |  2172 |  20.4624 |  4.0332 |  1.3423
11 |  800 |  $\\{0,\pm 1\\}$ |  1600 |  564 |  546 |  550 |  562 |  (8.46337) |  (0.4107) |  0.4775
14 |  800 |  $\\{0,1\\}$ |  4694 |  3064 |  3042 |  2993 |  3049 |  11.9894 |  (1.8703) |  0.0786
18 |  800 |  $\\{0,\pm 1\\}$ |  4694 |  992 |  983 |  942 |  986 |  20.4379 |  (8.418) |  1.52
22 |  2000 |  $\\{0,1\\}$ |  19990 |  13359 |  13340 |  13279 |  13352 |  2169.69 |  14.048 |  0.9654
27 |  2000 |  $\\{0,\pm 1\\}$ |  19990 |  3341 |  3316 |  3271 |  3339 |  2760.23 |  (23.568) |  27.768
32 |  2000 |  $\\{0,\pm 1\\}$ |  4000 |  1410 |  1364 |  1364 |  1398 |  (2909.38) |  (7.1411) |  32.101
35 |  2000 |  $\\{0,1\\}$ |  11778 |  7684 |  7656 |  7448 |  7634 |  1644.46 |  (14.806) |  1.2013

The computational performance of the SNN-CIM algorithm on GSET is also shown
in Table.1, which is compared with SA and CIM. We note that the algorithm of
SNN-CIM has an unparalleled advantage in computing time to find suboptimal but
adequate solutions for the application, both in the time and accuracy of the
solutions. Based on the above analysis, our proposed SNN-CIM has excellent
performance on combinatorial optimization problems.

##  III Conclusion

In addition, the model can be implemented by adding dissipative pulses to the
existing CIM based on fiber ring cavity, which is expected to improve its
computational performance and solve the problem of freezing in the excited
state during the calculation process. By using integrated photonics direct
coupling and high repetition frequency lasers, it is expected to increase the
pulse firing rate and the number of neurons, improve computing power, and
provide a better simulation platform for spiking neural networks. Finally, the
use of the signmoid function for the coupled term response in this model,
which is widely used in neural networks, creates an interesting connection
between CIM and neuromorphic computing, with the potential to excite the
research in the perceptron and pattern recognition of fiber-rings cavity based
CIM.

##  IV Acknowledgement

This research was funded by the National Natural Science Foundation of China
(Grants Nos. 62131002, 62071448), and the Fundamental Research Funds for the
Central Universities (BNU).

##  V Data availability

The GSET instances used for the tests is available at
https://web.stanford.edu/ yyye/yyye/Gset/. The SA algorithm temperature
adjustment is based on a logarithmic function, $T=T_{0}log_{2}(1+t_{0}/t)$.
All programs are written by matlab and run on intel i7-10700F 2.9GHz.

## References

  * Kitchen  _et al._ (2004) D. B. Kitchen, H. Decornez, J. R. Furr, and J. Bajorath, Nature reviews Drug discovery 3, 935 (2004).
  * Soler-Dominguez  _et al._ (2017) A. Soler-Dominguez, A. A. Juan, and R. Kizys, ACM Computing Surveys (CSUR) 50, 1 (2017).
  * Barahona  _et al._ (1988) F. Barahona, M. Grötschel, M. Jünger, and G. Reinelt, Operations Research 36, 493 (1988).
  * Arora and Barak (2009) S. Arora and B. Barak,  _Computational complexity: a modern approach_ (Cambridge University Press, 2009).
  * Kadowaki and Nishimori (1998) T. Kadowaki and H. Nishimori, Physical Review E 58, 5355 (1998).
  * Das and Chakrabarti (2008) A. Das and B. K. Chakrabarti, Reviews of Modern Physics 80, 1061 (2008).
  * Wang  _et al._ (2013) Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, Physical Review A 88, 063853 (2013).
  * Sutton  _et al._ (2017) B. Sutton, K. Y. Camsari, B. Behin-Aein, and S. Datta, Scientific reports 7, 1 (2017).
  * Saccone  _et al._ (2022) M. Saccone, F. Caravelli, K. Hofhuis, S. Parchenko, Y. A. Birkhölzer, S. Dhuey, A. Kleibert, S. Van Dijken, C. Nisoli, and A. Farhan, Nature Physics , 1 (2022).
  * Hamerly  _et al._ (2019) R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba, T. Honjo, _et al._ , Science advances 5, eaau0823 (2019).
  * Marandi  _et al._ (2014) A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Yamamoto, Nature Photonics 8, 937 (2014).
  * Inagaki  _et al._ (2016a) T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Yamamoto, and H. Takesue, Nature Photonics 10, 415 (2016a).
  * McMahon  _et al._ (2016) P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, _et al._ , Science 354, 614 (2016).
  * Inagaki  _et al._ (2016b) T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, _et al._ , Science 354, 603 (2016b).
  * Honjo  _et al._ (2021) T. Honjo, T. Sonobe, K. Inaba, T. Inagaki, T. Ikuta, Y. Yamada, T. Kazama, K. Enbutsu, T. Umeki, R. Kasahara, _et al._ , Science advances 7, eabh0952 (2021).
  * Fabian  _et al._ (2019) Fabian, Bhm, Guy, Verschaffelt, Van, der, and Sande, Nature communications 10, 3538 (2019).
  * Cen  _et al._ (2020) Q. Cen, T. Hao, H. Ding, S. Guan, Z. Qin, K. Xu, Y. Dai, and M. Li, arXiv preprint arXiv:2011.00064 (2020).
  * Pierangeli  _et al._ (2019) D. Pierangeli, G. Marcucci, and C. Conti, Physical Review Letters 122, 213902 (2019).
  * Fang  _et al._ (2021) Y. Fang, J. Huang, and Z. Ruan, Physical Review Letters 127, 043902 (2021).
  * Sun  _et al._ (2022) W. Sun, W. Zhang, Y. Liu, Q. Liu, and Z. He, Optics Letters 47, 1498 (2022).
  * Okawachi  _et al._ (2020) Y. Okawachi, M. Yu, J. K. Jang, X. Ji, Y. Zhao, B. Y. Kim, M. Lipson, and A. L. Gaeta, Nature communications 11, 1 (2020).
  * Leleu  _et al._ (2017) T. Leleu, Y. Yamamoto, S. Utsunomiya, and K. Aihara, Physical Review E 95, 022118 (2017).
  * Böhm  _et al._ (2021) F. Böhm, T. V. Vaerenbergh, G. Verschaffelt, and G. Van der Sande, Communications Physics 4, 1 (2021).
  * Reifenstein  _et al._ (2021) S. Reifenstein, S. Kako, F. Khoyratee, T. Leleu, and Y. Yamamoto, Advanced Quantum Technologies 4, 2100077 (2021).
  * Leleu  _et al._ (2019) T. Leleu, Y. Yamamoto, P. L. McMahon, and K. Aihara, Physical Review Letters 122, 040607 (2019).
  * Kako  _et al._ (2020) S. Kako, T. Leleu, Y. Inui, F. Khoyratee, S. Reifenstein, and Y. Yamamoto, Advanced Quantum Technologies 3, 2000045 (2020).
  * Leleu  _et al._ (2021) T. Leleu, F. Khoyratee, T. Levi, R. Hamerly, T. Kohno, and K. Aihara, Communications Physics 4, 1 (2021).
  * Inagaki  _et al._ (2021) T. Inagaki, K. Inaba, T. Leleu, T. Honjo, T. Ikuta, K. Enbutsu, T. Umeki, R. Kasahara, K. Aihara, and H. Takesue, Nature communications 12, 1 (2021).
  * Maass (1997) W. Maass, Neural networks 10, 1659 (1997).
  * Brette  _et al._ (2007) R. Brette, M. Rudolph, T. Carnevale, M. Hines, D. Beeman, J. M. Bower, M. Diesmann, A. Morrison, P. H. Goodman, F. C. Harris, _et al._ , Journal of computational neuroscience 23, 349 (2007).
  * Izhikevich (2003) E. M. Izhikevich, IEEE Transactions on neural networks 14, 1569 (2003).
  * Mnih  _et al._ (2015) V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, _et al._ , Nature 518, 529 (2015).
  * Yu  _et al._ (2012) Q. Yu, K. C. Tan, and H. Tang, in  _The 2012 International Joint Conference on Neural Networks (IJCNN)_ (IEEE, 2012) pp. 1–7.
  * Ma and Hao (2017) F. Ma and J.-K. Hao, Annals of Operations Research 248, 365 (2017).