#  A novel Three-step Network-based Ecosystem Modelling Framework

Laurin Steidle1, Inga Hense
_<EMAIL_ADDRESS>
_ Universitaet Hamburg, Grosse Elbstrasse 133, Germany_

(August 27, 2024)

###### Abstract

NEMF is a novel network-based ecosystem modelling framework. It is a flexible
and easy-to-use tool for modelling ecosystems with low- to intermediate
complexity. It is designed around the idea of visualizing an ecosystem through
a network that implicitly defines a set of differential equations. These
equations are then solved numerically and can also be inverse modelled. The
framework offers the functionality to handle non-equilibrium, non-linear
interactions.

## Introduction

We present a novel ecosystem modelling framework that is able to represent
non-linear interactions and offers inverse modelling capabilities. To our
knowledge it is the first framework that combines these features while still
being easy to use.

Research or management questions in ecosystem dynamics are often tackled
either with complex, specifically tailored ecosystem models or by deducing
trends from more general models that are well understood but do not represent
the studied system specifically.

However, both modelling approaches are limited. The first approach is highly
time consuming. The second inherently contains epistemic errors because the
more general model does not necessarily represent the system of interest
adequately.

Furthermore, when designing a new ecosystem model one needs to know its
interaction sufficiently well to estimate the parameters (i.e. mortality
rates) needed to represent them accurately. However, this is rarely the case.
Therefore, it is crucial to fit the unknown or unprecisely known parameters to
reproduce the dynamics of the observed ecosystem. This process is generally
referred to as inverse modelling or sometimes simply as optimisation or
fitting. In contrast, forward modelling is the process of using a model to
deduce predictions of a systems behavior. Solving the inverse modelling
problem is a time consuming challenge that comes on top of all the other tasks
required when modelling ecosystems.

Frameworks have been developed to simplify both forward and inverse modelling
problems. These provide tools for designing, running, fitting, and analyzing a
wide range of models for different applications. Famous examples in the marine
ecosystem community are FABM [1], carlibrar[2], and gadget [3] (We provide a
detailed overview of the presently available marine ecosystem frameworks
here.[4]) However, they all have major short comings.

  1. 1.

Some frameworks use linear differential equations [5, 6, 7, 8, 9] to describe
the interaction between system compartments. However, relationships within
ecosystems are highly nonlinear and thus this approach is too simplistic to
describe dynamic biological processes adequately.

  2. 2.

The frameworks with non-linear interactions are very complex, [1, 10] due to
the functionality they offer and their sheer size. This makes them too time
consuming to effectively apply them to less complex systems.

  3. 3.

The frameworks with non-linear differential equations and manageable
complexity are all designed for very-specific narrow applications [11, 3].
This makes them inapplicable for a all other ecosystem modelling tasks.

We propose a novel ecosystem modelling framework to overcome these
shortcomings. It is designed to offer an easy to use method for modelling
ecosystems with low- to intermediate complexity. The framework offers the
functionality to handle non-equilibrium, non-linear interactions. Because of
the broad research field of ecology and its research challenges, we designed
the framework in python with flexibility in mind. It provides a set of
functions for the typical ecosystem interactions but it is also possible for
the user to define their own interaction functions without needing to change
any of the framework code. The framework also offers methods that fit any
parameter of the model to represent the studied system.

For simplicity, the current version of the framework is limited to non-
spatially resolved models (box-models).

## Framework Overview

Due to the high diversity and complexity in ecosystems, there are typically no
established equations to describe them. Therefore, conceptual models are a
prerequisite for any ecosystem model development. Thus, typically the first
step when modelling an ecosystem is to create a schematic diagram to visualize
the ecosystem compartments and their interactions. The vertices in such a
network represent compartments in the system while the links also known as
edges represent their interactions. An exemplary network of an NPZD model
which is a simple marine ecosystem model is shown in fig. 1. (The so called
NPZD model, short for nutrient-phytoplankton-zooplankton-detritus model. See
e.g. [12] for a common representation.)

Figure 1: Network visualizing the interactions in an NPDZ model. Such a
network implicitly defines a a set of differential equations 1

Through such a type of network one implicitly defines a set of coupled
differential equations. In the case of the presented NPZD model in fig. 1 they
are given by

$\displaystyle\frac{\partial N}{\partial t}$
$\displaystyle=-J(N)P+\gamma_{m}D$ (1) $\displaystyle\frac{\partial
P}{\partial t}$ $\displaystyle=(J(N)-\Phi_{P})P-G(\epsilon,g,P)Z$
$\displaystyle\frac{\partial Z}{\partial t}$
$\displaystyle=(G(\epsilon,g,P)-\Phi_{Z})Z$ $\displaystyle\frac{\partial
D}{\partial t}$ $\displaystyle=\Phi_{P}P+\Phi_{Z}Z-\gamma_{m}D$

where the functions $J$,$G$,$\phi$,$\gamma$ represent the interactions, namely
growth, grazing, mortality, and remineralisation. The framework provides a
standardized way of defining the conceptual model through a simple human-
readable configuration file.

A major benefit of this approach is that all interactions can be easily
visualized and errors identified. This is particularly helpful when the
systems become more complex and therefore smaller errors might not become
immediately obvious.

In the second step, the differential equations are solved numerically. We
offer several numerical solvers for the time integration as there exists none
that suits all problems. Additionally, users may also define their time
integration methods.

Any parameter used in the forward model can also be inverse modelled. This may
include certain compartment values (i.e. amount of individuals) or interaction
parameters (i.e. mortality rates). To fit the model to a desired behaviour
(i.e. measurement data) one only needs to define the upper and lower bound for
the particular parameter in addition to providing the data describing the
desired behavior. The inverse problem is then expressed as a non-linear
programming problem with fixed constraints. It is solved by a set of different
solvers according to the user’s requirement. The conceptually most simple is
the so called gradient-descent approach discussed in the next section. A
description on how the framework is used including an example is found in the
next section while a more detailed one and further examples can be found in
the package Documentation [13].

## Description and Example

This section will use the previously introduced example of an NPZD model to
present the usage and technical background of the framework in detail.

### Forward modelling

####

The network shown in fig. 1 contains four vertices representing the
compartments. Each vertex is associated with a single value of a shared
quantity, often referred to as currency. The choice of currency varies
depending on the model. In our example we will use carbon mass in units of kg
as our currency.

Each interaction represents a process (e.g. grazing) that is modeled as a flow
from one compartment into another. There exist different types of flows. For
example, mortality based flows depend on one compartment exclusively, while
predator-prey interactions depend on both organism groups. To account for
this, each flow is labelled, representing the specific kind of interaction.

Furthermore, each vertex can be the origin and the destination of several
flows. The framework also allows that there might be multiple flows between
one set of origin destination vertices. Such an network is referred to as a
directed-labelled-multigraph. A complete example of a directed-labelled-
multigraph is shown in fig. 2.

Figure 2: A complete example of a directed multi-graph. a) shows a single
directed interaction. b) shows multiple directed interactions. c) shows a
bidirectional interaction.

We assume that the currency is conserved in all interactions. Hence, all the
flow leaving one vertex is always ending up in another. None is lost in the
transfer. Initially, this might seem like a limiting constrain. However, if
there is a desire to model an “inefficient” interaction where parts of the
flow are lost, an artificial loss compartment can be introduced.

####

To enable the framework to construct such kind of graph, the user only needs
to provide a list of flows, describing origin and destination, its type and
parameters used in the description of the flow. This set of information is
expressed in a standardized yaml format [14]. The interaction types must
either correspond to the name of an interaction function offered in the
frameworks interaction library or to a user written function.

The option to easily add user written alternatives to the present functions is
one of the core features of this framework as it offers a great deal of
flexibility which is not present in any other similar ecosystem framework.

####

A set of differential equations as shown in eq. 1 is constructed based on the
network structure read from the yaml file. The type of differential equations
that are constructed are non-linear first-order coupled ordinary differential
equations (ODE). These equations describe the dynamics of the system
deterministically for a given set of initial conditions. An exemplary
numerical solution of eq. 1 is presented in fig. 3

Figure 3: Exemplary solution of the differential equations shown in equation 1
as visualized in fig. 2. Each color represents one compartments population or
mass over time.

####

After the model has been designed and its numerical solution is found, the
question arises if this model actually represents the system that it is
designed to describe. Only rarely the first guess of the model parameters
define the model in such a way that it mimics the behaviour of the studied
system. A simple example of a system’s behaviour might be that the studied
system is in a steady-state. While the initial guess of the model might
perform reasonably well and also reaches a steady-state after a certain build-
up phase, it might reach a different one than the studied system.

For less complex models, certain parameters, often referred to as fitting-
parameters, can be adjusted by hand until the model mimics the studied system.
However, for more complex models this becomes unfeasible. Therefore, this so
called inverse modelling process needs to be automated.

####

When fitting the model by hand, its“quality” is implicitly judged. To automate
this process, an objective quality measure needs to be formalized. This
objective function is often also called loss- or cost-function. For our
framework we define the cost-function $L$ as the euclidean norm, which is a
quadratic measure of distance between the reference data provided for all
compartments i, and time t and the model output at these points.

$L(\alpha_{j})=||\bar{x}^{i}_{t}-\tilde{x}^{i}_{t}(\alpha_{j})||=\sqrt{\sum_{i,t}(\bar{x}^{i}_{t}-\tilde{x}^{i}_{t}(\alpha_{j}))^{2}}$
(2)

where $\bar{x}$ represents the reference value, and $\tilde{x}$ the model
value. $\tilde{x}(\alpha_{j})$ is a function of the model parameters
$\alpha_{j}$ which are tuned during the optimization process.

In a simple example there might only be two free parameters that are varied
during the optimisation process. In that case, each parameter would span one
dimension, while the cost at each points spans a third dimension as shown in
fig 4.

Figure 4: Exemplary cost-function field. Each point on this field represents
the measure of quality for a certain configuration of the model. The lower the
cost the better the model. Hence, it is the goal of the fitting process to
find the model configuration at the global minimum.

####

With the help of this cost-function field, which measures the quality of each
possible model configuration, a better model can be found through
minimization. Yet the key problem is, that the cost-function is not globally
known. Therefore, we have to explore the ”landscape” like a hiker in foggy
mountains. Each individual model configuration for a set of parameters has to
be evaluated to measure its performance. Note, that in general we can choose
any amount of free parameters creating a multi-dimensional cost-function
field.

This class of problem, where a continuous non-linear cost function is
minimized, is known as non-linear programming problem. There are several
different methods to tackle this problem.

####

One of the simplest method is a gradient based approach called ”gradient
descend”. [15] The key idea here is that as long as we always walk downhill we
will end up at the bottom of the valley eventually. To ensure that we always
walk downhill, we calculate the gradient at a given point by evaluating models
with slightly perturbed parameters. In other words, we take a step in each
direction and choose the direction that goes downhill the steepest. In its
simplest implementation, this concept is described by the following equation

$\alpha^{(n+1)}_{j}=\alpha^{(n)}_{j}-\lambda\nabla J(\alpha^{(n)}_{j})$ (3)

where $\alpha^{(n)}_{j}$ and $\alpha^{(n+1)}_{j}$ represents the current and
the ”improved” guess of model parameters. $\lambda$ describes the step size
that we take.

Typically one does not end up at the ideal model after the first step. Hence,
this process is iterated several times until either convergence is reached,
meaning that the model doesn’t improve anymore or a satisfying model is found.

A major problem in these optimization schemes is that generally, one can not
distinguish if the solution of the optimization process is a local or global
minimum. Meaning, that there might be a deeper valley on the other side of the
mountain that is unknown to us. (See fig. 4).

Unfortunately, there exists no guarantee that we find the deepest valley. The
only option we have to overcome this problem is by using different model
configurations as a starting point of the optimization process, and test if a
better solution can be found. This is known as a monte-carlo approach.

####

While the inverse problem is a challenging task from a mathematical
perspective, it is automated in our framework. Nevertheless, the user can
select several options for this process and also use self-written optimization
routines. The user only needs to provide two things:

  1. 1.

the data that represents the system that shall be represented

  2. 2.

to mark the parameters that are allowed to be varied by providing their upper
and lower bounds

####

After the framework has been configured and run, there are two possible kinds
of output depending on the configuration. In the case where no optimization
was applied, the framework returns the time evolution of the model. In the
second case, when the model was fitted to data, the time integration of each
iteration step, its associated parameters which configure the model, and costs
are returned. In both cases, a plot is automatically generated to visualize
the output. For the exemplary case, where we fit the NPZD model, the plot is
presented in fig. 5.

Figure 5: Exemplary output of a typical inverse modelling run. The lower
figure shows the fitted model with dashed lines representing the reference
data. The upper figure shows the model parameter configurations tested during
the fitting process.

The output can also be directly passed on for further processing, e.g. a
stability analysis.

## Discussion

The combination of the three steps:

  * •

conceptual design through a network

  * •

numerical integration of the underlying non-linear differential equations

  * •

inverse modelling

is a novel approach to ecosystem modelling which has not been used in any
ecosystem modelling framework that is not tailored to a very specific
application, i.e. fishery management. In the following we will discuss the
benefits and drawbacks of this method.

####

As discussed above, a major benefit of this framework is time that is saved
when using it compared to writing a model yourself. The configuration of the
framework for a large model, including its defining fitting constraints, can
be done in mere minutes. Through that the time that is needed from initially
drawing up a model, to the point where the results can be analyzed is reduced
to less then an hour. Writing a comparable implementation from scratch is
expected to take several days to a week. Therefore, this framework is ideally
suited for testing different implementations of nonlinear model interactions.

Due to its open design it is well suitable for a wide range of applications,
including management activities. This is not limited to fishery management as
many other frameworks but can also be applied in the area of e.g. blue carbon
management.

Because of its simple design is also well suitable for teaching. It can be
used with limited understanding of the inner workings. Hence, this allows it
to be used without in depth programming knowledge or understanding of the
numerics applied and therefore assisting in leading students to the field of
ecosystem modelling.

####

To this date, many ecosystem models and modelling frameworks used
approximated, linearized interactions to describe the studied ecosystem [5, 6,
7, 8, 9]. This is especially common for ecosystem-management applications. The
advantage of our approach is that it drastically simplifies the models and
therefore makes them much easier to use and interpret the results.

However, many crucial interactions such as predator-prey interactions are non-
linear, as already shown in the famous Lotka-Volterra equations. Hence,
linearized models are crude approximations that are valid if, and only if,
used to study systems with only small changes in the compartment quantities.
Therefore, a non-linear description of a system is necessary if its full
dynamics are studied.

####

The NEMF framework offers both forward and inverse modelling methods. However,
an inverse modelling framework does not necessarily need to offer an option to
solve the forward model. It can be designed in such a way, that the user can
provide an external forward model, for example, the calibrar[2] framework.
However, in that case the user needs to handle the communication between the
forward model that passes the cost function for a certain model configuration
to the inverse modelling framework. Therefore, they are required to write the
parser code themselves. This approach is sometimes referred to as _blackbox_
or _out-of-box_ approach.

In our framework, the forward modelling is included. This creates some
restriction concerning the model because the user is not free to provide their
own model anymore. However, it offers the large benefit that the user does not
need to write any code and does not need to understand how the inverse
framework is implemented, because they do not need to write any parsers. We
refer to this as _inside-the-box_ approach (See fig. 6).

Figure 6: Schematic representation of an inverse modelling framework. If the
upper forward modelling has to be provided by the user we call if a blackbox
approach while if it is included inf the framework we call it an inside-the-
box approach. NEMF is the latter.

####

While python offers many benefits, both for developing as well as for using
this framework, its computational speed is significantly slower compared to
e.g. C or FORTRAN. However, because we limit ourself to box-models this is not
an issue as computational times even for large single model runs is quasi-
instantaneous. Typically, full optimization runs with NEMF are performed in
time scales of seconds to several minutes for larger models.

The framework is designed for deterministic models exclusively. Therefore,
stochastic processes are not supported. However, we expect it to be possible
for the user to implement stochastic processes if needed. Yet, during the
development we do not test stochastic functions. Hence, we do not guarantee
that they work as expected.

####

While we only discussed the framework in the context of ecosystem modelling
applications, we want to mention that any kind of ODE system can be modelled
in the same way. We illustrated this in the framework documentation with the
example of an enzymatic reaction and coupled oscillated masses.

## Software Availability

The framework is published as open-source software under the BSD-3-Clause
License. The source code can be found on github.com/465b/nemf It can be
acquired as a python package and listed in the Python Package Index (PyPI)
under the acronym NEMF. For more information, see the official docs at
nemf.readthedocs.io.

## Acknowledgement

We want to thank Maike Scheffold for inspiring this work and offering her
feedback during the development. We also want to thank Jan Bydžovský and Klaus
Györgyfalvay for all the discussions that helped further this project.

## References

  * [1] Jorn Bruggeman and Karsten Bolding.  A general framework for aquatic biogeochemical models.  Environmental Modelling and Software, 61:249–265, 2014.
  * [2] Ramos Oliveros, Ricardo, Yunne, and Jai Shin.  calibrar: an R package for the calibration of complex ecological models, 2016.
  * [3] James Begley and D. Howell.  An overview of Gadget, the globally applicable area-disaggregated general ecosystem toolbox.  Ices Cm, pages 1–15, 2004.
  * [4] Overview of present marine ecosystem modelling frameworks.  https://465b.github.io/frameworks-overview/marine-ecosystem-modelling-frameworks.
  * [5] J. Melbourne-Thomas, S. Wotherspoon, B. Raymond, and A. Constable.  Comprehensive evaluation of model uncertainty in qualitative network analyses.  Ecological Monographs, 82(4):505–519, 2012.
  * [6] Stuart R. Borrett.  Throughflow centrality is a global indicator of the functional importance of species in ecosystems.  Ecological Indicators, 32:182–196, 2013.
  * [7] Stuart R. Borrett and Matthew K. Lau.  enaR: An R package for Ecosystem Network Analysis.  Methods in Ecology and Evolution, 5(11):1206–1213, 2014.
  * [8] Rowan Trebilco, Jess Melbourne-Thomas, and Andrew John Constable.  The policy relevance of Southern Ocean food web structure: Implications of food web change for fisheries, conservation and carbon sequestration.  Marine Policy, (January):103832, 2020.
  * [9] N. Niquil, B. Saint-Béat, G. A. Johnson, K. Soetaert, D. van Oevelen, C. Bacher, and A. F. Vézina.  Inverse Modeling in Modern Ecology and Application to Coastal Ecosystems, volume 9.  Elsevier Inc., 2012.
  * [10] Asta Audzijonyte, Heidi Pethybridge, Javier Porobic, Rebecca Gorton, Isaac Kaplan, and Elizabeth A. Fulton.  Atlantis: A spatially explicit end-to-end marine ecosystem model with dynamically integrated physics, ecology and socio-economic modules.  Methods in Ecology and Evolution, 10(10):1814–1819, 2019.
  * [11] Yunne Jai Shin and Philippe Cury.  Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing.  Canadian Journal of Fisheries and Aquatic Sciences, 61(3):414–431, 2004.
  * [12] Peter J.S. Franks.  Coupled physical‐biological models in oceanography.  Reviews of Geophysics, 33(2 S):1177–1187, 1995.
  * [13] NEMF documentation.  https://nemf.readthedocs.io.
  * [14] YAML documentation.  https://yaml.org/spec/1.2/spec.html.
  * [15] Haskell B. Curry.  The method of steepest descent for non-linear minimization problems.  Quarterly of Applied Mathematics, 2(3):258–261, 1944.


# When costly migration helps to improve cooperation

Hsuan-Wei Lee  Colin Cleveland  Institute of Sociology, Academia Sinica,
Taiwan  Attila Szolnoki  Institute of Technical Physics and Materials Science,
Centre for Energy Research, P.O. Box 49, H-1525 Budapest, Hungary
<EMAIL_ADDRESS>
###### Abstract

Motion is a typical reactions among animals and humans trying to reach better
conditions in a changing world. This aspect has been studied intensively in
social dilemmas where competing players’ individual and collective interests
are in conflict. Starting from the traditional public goods game model, where
players are locally fixed and unconditional cooperators or defectors are
present, we introduce two additional strategies through which agents can
change their positions of dependence on the local cooperation level. More
importantly, these so-called sophisticated players should bear an extra cost
to maintain their permanent capacity to evaluate their neighborhood and react
accordingly. Hence, four strategies compete, and the most successful one can
be imitated by its neighbors. Crucially, the introduction of costly movement
has a highly biased consequence on the competing main strategies. In the
majority of parameter space, it is harmful to defectors and provides a
significantly higher cooperation level when the population is rare. At an
intermediate population density, which would be otherwise optimal for a system
of immobile players, the presence of mobile actors could be detrimental if the
interaction pattern changes slightly, thereby blocking the optimal percolation
of information flow. In this parameter space, sophisticated cooperators can
also show the co-called Moor effect by first avoiding the harmful vicinity of
defectors; they subsequentially transform into an immobile cooperator state.
Hence, paradoxically, the additional cost of movement could be advantageous to
reach a higher general income, especially for a rare population when subgroups
would be isolated otherwise.

What should we do if we want to avoid being exploited by others or to realize
a more supportive environment? A simple answer is to leave our location and
search for a better place. This idea was studied intensively in evolutionary
game theoretical models where cooperation and defection compete. The problem
has been investigated from different angles, where an individual’s ambition to
leave a detrimental place or search for a more promising location could be the
main motivation. Interestingly, the majority of previous models assumed
homogeneous populations where everyone is ready to move, ignoring the extra
investment required to permanently monitor our neighborhood. In the present
research, we explore both aspects and assume that staying or thinking about
leaving could be competing options, with explicit consideration to the
additional cost of the latter. In the resulting 4-strategy model, we provide
insight into the conditions under which costly migration could be beneficial
to attain higher general well-being across the entire society.

##  I Introduction

Wars, climate change, and decaying living conditions are the major factors
that spur migration, which not only hugely affects the global economy, but
also precipitates further change in our hectic world. In 2019, for instance,
at least 3.5 % of the world’s population comprised international migrants,
according to the UN McAuliffe and Khadria (2019). Individual mobility,
however, is not just a simple reaction to a specific living condition; it
could also be the source of a co-evolutionary mechanism of how collective
well-being evolves in the framework of human cooperation Perc and Szolnoki
(2010); Xu and Hui (2022); Ohdaira (2022); Pi  _et al._ (2022); Kang  _et al._
(2021); Tao  _et al._ (2021); Couto  _et al._ (2022); Liu and Wu (2022).
Therefore, when we sought to understand the evolution of cooperation in
spatially structured populations, considering mobile agents was among our core
ideas from the outset Kendal (1992); Vainstein  _et al._ (2007); Sicardi  _et
al._ (2009).

Previous models have focused on different aspects of individual motivations to
move. An involved actor’s reaction could be risk-driven Chen  _et al._ (2012)
or environment-driven Xiao  _et al._ (2022); Park  _et al._ (2013); Wang  _et
al._ (2022); Park (2019), but aspiring to a declared income or expectation
regarding a certain payoff level could also be a prime goal of migration Lin
_et al._ (2011); Wu  _et al._ (2012); Chen  _et al._ (2016a). The principal
question along this research avenue is how migration influences the evolution
of cooperation among players who are involved in a social dilemma where their
individual and collective interests are in conflict. More precisely, they can
collect a higher payoff individually if they defect, but if they all choose
this strategy, then they will earn a smaller income than in the scenario when
they cooperate. Interestingly, the conclusions are more subtle in terms of
summary in the form of a simple message Cong  _et al._ (2017); Park (2022); Li
_et al._ (2020); Chen  _et al._ (2016b); Xiao  _et al._ (2020); Dhakal  _et
al._ (2022). In some cases, the players’ movement is proved to be harmful,
while in other cases, it clearly becomes a mechanism for supporting
cooperation. Additionally, there are examples where the intensity of mobility
is a decisive factor in the evolutionary outcome Lütz  _et al._ (2021);
Cardinot  _et al._ (2019).

Nevertheless, previous model studies typically assumed that a homogeneous
population where every player has the capacity to move and improve their
individual payoff. This choice could be the result of a decision given that
moving could be a costly skill. More precisely, those players who are willing
to move should constantly watch and evaluate their neighborhood, as whether to
move will depend on the result of this analysis. Undoubtedly, these efforts
require an extra investment from the mentioned players, which should be
considered an extra cost. Perhaps it is worth mentioning that one previous
study examined costly migration Liu  _et al._ (2012), but the conclusions were
controversial. To solve the puzzle and provide a more comprehensive conclusion
in our research, we assume that the population could be heterogeneous. More
precisely, a player can decide whether she wants to invest extra efforts to
explore the neighborhood and possibly change residence, or whether,
alternatively, she prefers to save in terms of the payoff by not moving and
using unconditional cooperator or defector strategies. Premised on this, we
constructed a 4-strategy model where both approaches are available, allowing
us to explore the profitability of each attitude at specific parameter values,
namely, average population density and the strength of the social dilemma. In
the next section, we define our extended public goods game (PGG) and introduce
the key parameters. Sec. III highlights our main observations. Sec. IV
concludes with a summary of the results and a discussion about their
implications.

##  II Model

We start from a traditional spatial PGG, where players are distributed on a
square lattice with a linear size of $L$, and periodic boundary conditions are
applied Sigmund (2010); Perc  _et al._ (2013). Some of the lattice points are
occupied by a player while others are not; hence, altogether, $0<\rho\leq 1$
portion of the available $L\times L$ positions are filled. According to the
standard protocol, an actor, who could be an unconditional cooperator ($C$) or
an unconditional defector ($D$), plays a PGG with the nearest available
neighbors in which a cooperator player invests the amount of $c=1$ in the
common pool, while a defector player refrains from doing so. We then increase
the accumulated contributions by an $r$ enhancement factor and distribute the
resulting public goods equally amongst all group members regardless of their
strategies. Of course, players are also involved in the games their neighbors
organize.

We therefore extend the model by introducing two new strategies that behave
similarly to the original strategies in the game but entail an extra skill
related to moving. More precisely, it is well-known that a cooperation level
that is too low is detrimental to everyone. Therefore, the new group of
players, called sophisticated players, are capable of evaluating their local
neighborhood and leaving it if it is harmful. In particular, they will leave
their original positions if the average cooperation level is less than 50%
among the existing neighbors. Moreover, if there is no other actor in the
nearest neighborhood because all four neighboring lattice sites are empty, the
actor will leave the site, too. In this way, we introduce sophisticated
cooperator ($SC$) and sophisticated defector ($SD$) strategies. These players
always explore their environment, and if they are dissatisfied, they move to a
randomly chosen neighboring empty place and make a trial strategy change in
the new position. Importantly, sophisticated players bear an extra cost $f$
because they need to constantly watch and evaluate their neighborhood. One may
argue that the introduction of mobile players will definitely support
cooperation because cooperators may have a chance to escape from a situation
where they are being exploited. On the other hand, mobile defectors also have
the opportunity to search for remote players who could be new prey given their
selfish attitude. Hence, the question is far from trivial. Interestingly,
however, as we will show, this extra cost has a biased consequence that
impacts the competing strategies.

The strategy change is based on the widely used imitation process, where the
probability of a player $x$ adopting the strategy of a neighboring player $y$
is based on the payoff difference between these competitors:

$W(s_{x},s_{y})=\frac{1}{1+\exp[(\Pi_{x}-\Pi_{y})/K)]},$ (1)

where parameter $K$ quantifies the uncertainty of strategy imitation. To
attain results that are comparable with those previously obtained in the
traditional model, we fixed $K=0.1$.

In sum, at the beginning of the Monte Carlo (MC) simulation, we first randomly
distribute $N=\rho\times L\times L$ players on a square lattice where $\rho$,
as a key parameter, determines the average population density. The players
randomly choose a strategy from the available set $(C,D,SC,SD)$. During an
elementary step, we randomly choose a player $i$. If the strategy that player
$i$ selects is $SC$ or $SD$, which means we have a sophisticated player, then
player $i$ explores her neighborhood and determines the average cooperation
level among the neighbors. If it is below $1/2$, then the player moves to a
randomly chosen empty neighboring place. In the next step, player $i$, either
in the original position or in the new place, selects a neighboring player $j$
who could be a potential model player for strategy adoption. After we
calculate the payoff values for both involved players according to the PGG
protocol, it is reduced by cost $f$ for a sophisticated player. Next, player
$i$ imitates player $j$ with the probability defined by Eq. 1. Evidently,
traditional players are only involved in the strategy updating step at their
original positions. A full MC step is established if we repeat the above-
described procedure $N$ times. The relaxation time to reach the stationary
state depends not only on the applied system size but also on the population
density. In each case, we monitored the evolution of strategies and measured
the proportions of strategies when they reached the stationary values. The
typical relaxation time was about 10,000 MC steps. In this research, we
typically applied a linear system size of $L=400$, but we also checked other
sizes to exclude finite-size effects. To obtain reliable values for the
cooperation level or for the proportions of strategies, we have typically
averaged the results of 10-50 independent runs.

##  III Results

Let us begin by presenting our key observation about the consequence of mobile
players’ presence in addition to traditional immobile actors. This can be done
effectively if we compare the general cooperation level for the original
2-strategy model with that of the extended 4-strategy case defined in the
previous section. As we noted, our main control parameters are the $r$ synergy
factor of PGG, which characterizes the strength of the dilemma, and $\rho$
general density, which indicates the population’s degree of rarity or
crowdedness.

A comparison of the results obtained in the mentioned cases is presented in
Fig. 1. The upper panel shows the traditional model, where we can see that
there is an optimal density that provides the highest cooperation level,
especially at modest $r$ values. This result confirms a previous finding that
at a fixed dilemma strength, we can reach a higher cooperation level if the
population is less crowded Wang  _et al._ (2012a). More precisely, the highest
cooperation level can be reached in the vicinity of the percolation threshold,
which provides a successful combination of external conditions: The population
is sufficiently rare to avoid providing easy prey for defectors, but on the
other hand, information about cooperators’ collective success can spread
throughout the whole system Wang  _et al._ (2012b). The latter condition is
not fulfilled if the population is rare. In this case, it is almost a random
situation in which strategy wins on a small island and the local evolutionary
outcome depends largely on the initial distribution of strategies. Therefore
the final and practically frozen state is a mixture of pure $C$ and pure $D$
islands. This consequence can be seen nicely in the upper panel of Fig. 1,
where the cooperation level is around 0.5, almost independently of the dilemma
strength if $\rho$ is sufficiently small.

Figure 1: Cooperation level on the enhancement factor–density parameter plane
for two models. The upper panel shows the result of the 2-strategy model where
unconditional cooperators ($C$) and defectors ($D$) are present and players
have no moving capacity. The bottom panel shows the cooperation level for our
present 4-strategy model, which represents our extension of the original model
with the addition of two new strategies that entail cooperation ($SC$) or
defection ($SD$), and players are capable of leaving their positions if they
find them detrimental. This skill requires extra effort, in the form of an
extra fee, which is $f=0.08$ in the plotted case. Evidently, the sum of the
proportions of $C+SC$ fractions indicates the cooperation level in the latter
case.

This situation changes significantly if mobile strategies are also introduced.
In this case, a significantly higher cooperation level is also reached in a
rare population. One might say that the information spreads due to mobile
agents, but the explanation is not as simple because defectors can also move.
The deeper answer is based on the highly biased consequence of costly
movement, which is harmful to defectors but can offer cooperators extra help.
We will elaborate on this argument later. Here, it is worth noting that the
applied extra fee $f$ is sufficiently high to represent the enhanced
investment for sophisticated players, but it is not too large to make their
payoff values unfeasible.

If we carefully compare the mentioned panels, we can see that a better outcome
can be reached in the traditional 2-strategy model around the previously
mentioned optimal $\rho\approx 0.6$ density when $r$ is small; hence, the
dilemma strength is more critical. For a quantitative comparison, in Fig. 2,
we simultaneously present both models’ cooperation levels dependent on the
synergy factor at a fixed value of $\rho=0.65$. Evidently, in the 4-strategy
model, the cooperation level is calculated as the sum of the proportions of
the $C$ and $SC$ strategies. The figure shows that a significantly higher
cooperation level can be reached in the traditional model, where players
cannot move, and the difference only vanishes for higher $r$ values. This
effect can be understood if we consider how movements happen. If a player is
satisfied because the cooperation level is sufficiently high, then even mobile
players will maintain their positions and move only in the alternative case.

Figure 2: Cooperation level dependent on the enhancement factor at $\rho=0.65$
for the traditional 2-strategy model (red) and the extended 4-strategy model
(blue), with $f=0.08$, as indicated. Around the optimal density, which is
close to the percolation threshold, mobile actors’ presence could be harmful
when the synergy is low.

This so-called environment-driven movement can cause slight changes in the
interaction graph, resulting in more aggregated colonies. The aggregation of
players can be measured directly if we calculate the probability of finding
nearest-neighbor players around a player. In this way, two-point probability
$p(1,1)$ can indicate the likelihood of finding the nearest neighbor next to
our focal player. For example, if we consider the two-strategy model at the
mentioned density of $\rho=0.65$, then the supposed random distribution of
players results in an approximate value of $p(1,1)=0.42$ for the mentioned
pair probability. This value is significantly larger for the 4-strategy model,
where players may move and approach each other. More specifically, if we use
the same $\rho$ value at $r=3.4$ and $f=0.08$, then the mentioned probability
increases to 0.60. In other words, the population becomes more aggregated, but
it also means that the comprehensive percolation is in jeopardy because we are
close to the critical $\rho$ value. In this way, mobility may threaten one of
the conditions we previously identified as necessary to reach the highest
cooperation level. Of course, if the synergy factor is sufficiently high, then
we do not need proper percolation because the cooperator can also win locally;
hence, the detrimental consequence of the player’s movement becomes
irrelevant, and the cooperation level is restored to reach the value produced
by the simpler model with fixed players. To summarize our observations, there
is a small area on the parameter plane at intermediate $\rho$ values and
modest $r$ values when the traditional model provides a better outcome, while
for the majority of parameter pairs, the 4-strategy model undoubtedly offers a
higher cooperation level.

In the following paragraphs, we provide deeper insight into the consequences
of adopting new strategies. In Fig. 3, we present a phase diagram on the
$r-\rho$ parameter plane, where the same extra cost value for sophisticated
players is used as in Fig. 1. First, we would like to stress that only
traditional defectors survive if the synergy factor is sufficiently low,
independent of the average density $\rho$. Evidently, if the social dilemma is
strong, then cooperators cannot gain the necessary income to compensate for
defectors’ exploitation. Consequently, only those who adopted $D$ and $SD$
strategies survive at the intermediate evolutionary stage. In this case,
however, the extra cost of the $SD$ strategy forms a disadvantage comparing
with simple $D$ players, who will eventually dominate the whole system.

Interestingly, the advantage of the $D$ over the $SD$ strategy remains even if
the extra cost $f$ is zero. In the latter case, one might expect that the
capacity to move should provide some benefit to adopters of the $SD$ strategy,
but the opposite is true. The reason is as follows. The advantage of $SD$ to
move also makes accumulation harder; it is therefore easier for an $SD$ player
to meet a $D$ player in a new position. On the contrary, if we select a $D$
player, who cannot escape from nearby players for strategy update, then it is
more likely that the $D$ player will choose another $D$ neighbor as a model
player, in which case nothing happens, and the portion of $D$ players remains
unchanged. When we select an $SD$ player during the elementary process, then
she may move first, but she may sequentially encounter a $D$ player. Indeed,
there is no difference between their payoff values, but Eq. 1 indicates a 0.5
probability of adopting strategy $D$, which can decrease the portion of $SD$.
This asymmetry will gradually lead to a fully $D$ destination.

Figure 3: Phase diagram of the 4-strategy model on the $r-\rho$ plane obtained
for an additional cost $f=0.08$. If the synergy factor is sufficiently low,
then only fixed defectors survive, independent of the average population
concentration. Sophisticated defectors can only survive in a limited parameter
area if the dilemma is mild and the system is crowded. Sophisticated
cooperators, however, are viable in a large area, especially when mobility has
an enhanced role in a rare population.

Returning to the phase diagram, it highlights that the $SD$ strategy is only
viable in a very limited parameter area when the dilemma strength is modest
and the population is practically crowded. At these $r$ values, the
coexistence of cooperator and defector strategies is expected in a spatial
system Szolnoki  _et al._ (2009); therefore, $SD$ players can enjoy the
proximity of supporting cooperators that can generously cover their extra
fees. On the other hand, their skill of leaving a detrimental environment
makes them more competitive relative to the pure $D$ strategy. In other words,
our preliminary expectation about the advantage of mobile players is justified
in the mentioned case. Importantly, however, it is only valid for the $SD$
strategy, not the $SC$ strategy. Players who adopt the latter strategy simply
have no chance of escaping defectors because of the high general density,
making it evident that the extra cost puts them at a disadvantage. These
mechanisms explain that $C+SD$ forms a stable solution here.

When the system is less crowded, another aspect of the actor’s mobility can be
detected. In a rare population, cooperators may have a good chance of escaping
defectors; hence, $SC$ strategy becomes viable. Interestingly, intermediate
concentration does not provide an opportunity for the $SD$ strategy. Indeed,
these $SD$ adopters can escape from nearby players, but they can no longer
utilize cooperator strategies, because the latter type of player can form a
compact island, which gives them highly competitive payoff values. As we
stressed above, intermediate density was already a perfect condition in the
2-strategy model, and this conclusion remains intact in the extended model.
One may claim that $C$ and $SC$ are still unequal since the latter should bear
the extra cost $f$. However, they remain separate because $SC$ players can
form “satisfied” groups; hence, their coexistence is a mixture of the $C$ and
$SC$ domains.

Finally, if the population is very rare, then immobile defectors once again
have a chance at survival. In this case, they win their local battle, but they
are stuck afterward. Other domains invaded by $C$ or $SC$ players are also
satisfied; they remain still, and the information about more successful
strategies does not reach $D$ players. This explains why we can observe the
$D+C+SC$ solution here.

Before discussing these parameter regions, let us return to the case of
intermediate $\rho$ values; we can observe an interesting process if we
monitor the proportion of strategies during the evolution. Typical time
evolution is shown in Fig. 4, where we plotted how the proportions of
strategies change if we launch the evolution from a random initial state at
$r=4$, $\rho=0.65$, where sophisticated players bear an extra cost of
$f=0.08$. At the very beginning, for a short initial period, defector players
enjoy a random local neighborhood; hence, cooperators are decaying. The
survivors then form cooperative clusters that provide them with competitive
payoffs; hence, their numbers start growing. This “first down, later up”
dynamic is the typical trademark of network reciprocity reported earlier in
several cases Perc  _et al._ (2008); Szolnoki and Perc (2009); Wei  _et al._
(2021); Shen  _et al._ (2022); Liang  _et al._ (2022). Another consequence of
the emergence of cooperative clusters is the rapid extinction of the $SD$
strategy. Even if they can escape nearby players, during their migration, they
may encounter a phalanx of cooperators that can easily defeat them. Of course,
this effect is not found among pure $D$ players, who cannot move; hence, they
are not threatened by a more successful strategy. Accordingly, their
proportion declines at a significantly slower rate due to the strategy
invasion process. The portion of sophisticated cooperators, however, starts
declining again after reaching the maximum level. This phenomenon can be
explained as follows. First, it is advantageous for cooperators to be mobile
because they can avoid intensive exploitation at the hands of defectors. This
is why their minimal proportion in the first stage is a bit higher than the
minimum level for the $C$ strategy. Hence, the presence of the $SC$ strategy
helps to fight defection more effectively. However, later, when defectors are
eradicated, the vicinity of the simple $C$ strategy becomes dangerous to them
because the negative consequence of the additional cost $f$ lowers their
payoff. This process resembles the so-called Moor effect Szolnoki  _et al._
(2011). When sophisticated cooperators have performed their duties (and
weakened defectors), they may go (for the benefit of simple cooperators).

Figure 4: The early evolution of the proportions of strategies in the extended
model at $r=4$, $\rho=0.65$, and $f=0.08$. We can detect three different
stages in the evolutionary process, which are discussed thoroughly in the main
text. While the early growth of $C$ and $SC$ strategies after an initial
decline is a sign of how network reciprocity works, the subsequent decline of
the $SC$ strategy is a consequence of the so-called Moor effect. Finally,
though this is not shown here, all defector strategies die out, and only
cooperator strategies survive. Figure 5: Representative strategy distributions
obtained for a density of $\rho=0.65$ at $r=3$ (left column) and at $r=5$
(right column). Panels (a) and (b) show the 2-strategy model, while panels (c)
and (d) depict the 4-strategy case. The panels only show a $100\times 100$
section of a larger system. Dark (light) blue represents pure (sophisticated)
cooperators, while dark red shows pure defectors. Note that sophisticated
defectors cannot survive here.

In Fig. 5, we present some characteristic patterns in this parameter area to
summarize the differences between the traditional and the extended models. The
first row shows the stationary strategy distribution in the 2-strategy model
while the bottom row depicts the results of the 4-strategy case. Furthermore,
the patterns shown in the first column were obtained at $r=3$, which
represents harsh external conditions, while the second column illustrates
observations for a more friendly case at $r=5$. The difference between the two
is striking. While the players’ spatial distribution is random in panels (a)
and (b) due to the model’s similarity to the traditional 2-strategy model, the
players are more aggregated when individual mobility is allowed. This is
especially true for a small $r$, as shown in panel (c), because in this case,
nobody is satisfied in the initial stage; hence, sophisticated players move
actively. However, they are prone to easy capture by defectors because of the
low $r$ value. In the large $r$ case, the population is still more aggregated
than a random distribution; hence, the pattern is below the percolation
threshold, though the total number of players would dictate the percolation in
a uniformly distributed system. This explains why $SC$ players can survive:
because they can form isolated domains comprising pure $C$ players who would
otherwise beat them. For a similar reason, some isolated $D$ spots remain, but
they represent a negligible minority; hence, the system closely approaches a
full cooperator state.

Next, we discuss the low $\rho$ region when the population is rare and smaller
communities become partly separated. In this case, the introduction of costly
mobile strategies has a clear positive impact regarding improving general
cooperation. This can be seen in Fig. 6, where we compared the average
cooperation levels dependent on the synergy factor at $\rho=0.25$. First, we
briefly discuss the gentle increasing red slope of the 2-strategy model. As we
mentioned at the beginning of Sec. III, here, the key parameter that
characterizes the dilemma strength only has second-order importance. There are
small communities where the evolutionary outcome is sensitive to the initial
distribution of strategies. Therefore either cooperators or defectors can win
locally, which provides an average cooperation level close to the random
value. Of course, the value of the synergy factor plays a certain role; hence,
the cooperation level increases slightly as we increase $r$. However, the
$r$-dependence is unusually weak, which is the consequence of poor
communication between smaller communes. Of course, it also means that
cooperators may survive at very small $r$ values, but practically, it is an
artificial effect, and it is divorced from the social dilemma.

Figure 6: Cooperation level dependent on the enhancement factor at $\rho=0.25$
for the traditional 2-strategy model (red) and the extended 4-strategy model
(blue), with $f=0.08$, as indicated. In this parameter region, the presence of
sophisticated players has a clear positive impact on public cooperation.

This situation changes dramatically when we introduce mobile players. In this
case, a significantly higher cooperation level can be reached at a low $r$
value because mobile players can spread the good news about the benefits of
cooperation. Of course, $SD$ players can also move, but they are unprotected,
as we argued earlier. If we compare the blue line with the results shown in
Fig. 2, we see that full cooperation cannot be reached in a rare population,
even if we have mobile actors. This is because some tiny domains are occupied
by pure $D$ players. Even if they are unsatisfied, they cannot change their
status. Other domains formed by $C$ players are also immobile, or, if they are
$SC$ players, then they are satisfied in the company of nearby neighbors.
Therefore, they are not motivated to move; hence, the poor performance of
defector spots remains hidden, which noticeably lowers the general cooperation
level. Interestingly, this mechanism cannot work at an intermediate density
because defector domains cannot be completely isolated. It is also worth
noting that the above-mentioned Moor effect regarding the time evolution of
the $SC$ strategy is less visible at small densities, the explanation for
which is similar to that for the survival of fixed defectors. Specifically,
isolated $SC$ actors may reach a state of satisfaction in the vicinity of
similar players; hence, they do not move, and can keep their original strategy
because their additional cost no longer has any relevance.

Figure 7: Representative strategy distributions obtained for a density of
$\rho=0.1$ at $r=3$ (left column) and at $r=5$ (right column). Panels (a) and
(b) show the 2-strategy model, while panels (c) and (d) depict the 4-strategy
case. The panels only show a $100\times 100$ section of a larger system. The
color code is the same as in Fig. 5.

Next, we support our arguments by presenting some representative patterns
obtained at $\rho=0.1$. In Fig. 7 we follow the same logic as in Fig. 5. As
the top row highlights, we have small isolated homogeneous domains in the
final state, where either $C$ or $D$ wins. If we compare the columns, there is
no significant detectable difference between the configurations shown in
panels (a) and (b). Here, the $r$ value has no crucial role because the final
local outcome is highly sensitive to the initial distribution of strategies.

Turning to the 4-strategy case, $SD$ cannot survive; hence, we did not detect
players represented by light red. Compared to the top row, the aggregation of
players is obvious in panels (c) and (d) as a result of sophisticated players’
migration. Due to the small $r$ value, a small island should be homogeneous.
Here, either cooperation or defection should prevail locally. Notably, both
$C$ and $SC$ strategies can survive. Adopters of the latter may reach a state
of satisfaction in the vicinity of similar actors and therefore refrain from
any further movement. For the high $r$ case, the main difference is that
defectors cannot build large domains it is no longer as easy for them to
emerge victorious from the initial battle. This mechanism results in a
significantly higher general cooperation level. However, some isolated red
spots remain, which explains why we cannot realize a full cooperator state
here. The latter would require an effective information percolation throughout
the whole system.

Finally, we briefly discuss the role of the additional cost parameter $f$,
which is a crucial element of our model compared to the large number of other
models previously used to understand the impact of migration on the original
social dilemma. Evidently, if this cost is too high, then sophisticated
strategies become less viral, and we practically revert to the original
2-strategy model. Interestingly, however, the behavior of the $SD$ strategy is
almost independent of the value of $f$, and the strategy only survives in a
very limited parameter area when $\rho$ is sufficiently high and $r$ is in an
intermediate region. From this aspect, the presented behaviors obtained at
$f=0.08$ faithfully represent what we can expect of other moderate $f$ values.

##  IV Conclusion

It is perhaps unnecessary to emphasize the importance of exploring the
possible impacts of migration on the outcomes of social dilemmas because
individual mobility may generally result in subtle collective behavior among
system members Vicsek and Zafeiris (2012). Therefore, the last decade has seen
intensive research activity along this path Zheng  _et al._ (2022); Flores
_et al._ (2021); Quan  _et al._ (2021a); Deng and Zhang (2022); Lee  _et al._
(2022); Liu  _et al._ (2021); Lv and Song (2022); Quan  _et al._ (2021b);
Zheng  _et al._ (2021); Li  _et al._ (2021a); Fu  _et al._ (2021); Zhu  _et
al._ (2021); Yang and Zhang (2021); Li  _et al._ (2021b); Lee  _et al._
(2021); Quan  _et al._ (2020); Fu  _et al._ (2019). In this research, our
principal motivation was twofold. First, we wanted to ascertain the impact if
players have the liberty to decide whether they want to ponder the option of
moving. This aspect can be modeled using two significantly different families
of strategies. While the first one focuses on how to behave in the original
social dilemma, the other entails permanent observation of the neighborhood,
and individuals may move if their neighbors’ behavior is unsatisfactory. In
other words, it could also be the source of a dilemma regarding whether we
want to ponder with the option of moving. If the answer is yes, then we will
pay the price for extra information and for maintaining the capacity to change
our position. Alternatively, we can consider the cost of proper movement into
this parameter. Notably, the majority of previous work has ignored this
aspect.

It is also worth mentioning that we introduced the so-called sophisticated
players as strategy-neutral, because both cooperators and defectors may choose
this option. Therefore, its impact for the cooperation level is far from
trivial. Our main observation was that the introduction of an extra cost has a
highly biased effect on the competition among the original strategies. In
particular, the sophisticated cooperator strategy is viral and can help to
elevate the cooperation level in the majority of the available area within the
parameter space. This is especially true in a population where the average
density is low and information exchange would be difficult otherwise. There is
a very limited parameter space, with an intermediate concentration and a low
synergy factor, where sophisticated players’ presence is detrimental. In this
narrow parameter region, the fixed players would have formed an optimal
interaction pattern, which provides a high cooperation level. In the mentioned
case, the population is sufficiently rare; hence, defectors have little
opportunity to exploit neighboring cooperators. However, because we are beyond
the percolation threshold, the positive impact of cooperation can spread
throughout the whole system. Allowing some players to move damages the latter
condition. Consequently, players aggregate more intensively, and the
percolation is broken.

The asymmetric impact of the extra cost on strategies is also evidenced by the
fact that sophisticated defectors can only survive in a very limited parameter
area in coexistence with pure cooperators. In this section of the parameter
plane, the population is essentially crowded; hence, individual mobility only
has second-order importance. Furthermore, the applied synergy factor would
also allow the coexistence of pure strategies.

In sum, our model is an initial trial to investigate more realistic models
where all aspects of individual mobility are considered, allowing us to
accurately estimate the impact of this factor on fundamental social dilemmas.
Extensions to other evolutionary games and alternative interaction graphs
would increase our understanding of the original question.

## Acknowledgments

The research reported was supported by the Ministry of Science and Technology
of the Republic of China (Taiwan), under grant No. 109-2410-H-001-006-MY3.

## Author declarations

The authors declare no competing interests.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## References

  * McAuliffe and Khadria (2019) M. McAuliffe and B. Khadria, eds.,  _World Migration Report 2020_ (Int. Org. for Migration, Geneva, Switzerland, 2019).
  * Perc and Szolnoki (2010) M. Perc and A. Szolnoki, BioSystems 99, 109 (2010).
  * Xu and Hui (2022) C. Xu and P. M. Hui, Phys. Rev. E 105, 054309 (2022).
  * Ohdaira (2022) T. Ohdaira, Sci. Rep. 12, 6604 (2022).
  * Pi  _et al._ (2022) B. Pi, Z. Zeng, M. Feng, and J. Kurths, Chaos 32, 023117 (2022).
  * Kang  _et al._ (2021) H. Kang, X. Zhou, Y. Shen, X. Sun, and Q. Chen, Phys. Lett. A 417, 127678 (2021).
  * Tao  _et al._ (2021) Y. Tao, K. Hu, and L. Shi, EPL 135, 28001 (2021).
  * Couto  _et al._ (2022) M. C. Couto, S. Giaimo, and C. Hilbe, New J. Phys. 24, 063010 (2022).
  * Liu and Wu (2022) F. Liu and B. Wu, Appl. Math. Comput. 431, 127309 (2022).
  * Kendal (1992) J. K. Kendal, Evolution 46, 1492 (1992).
  * Vainstein  _et al._ (2007) M. H. Vainstein, A. T. C. Silva, and J. J. Arenzon, J. Theor. Biol. 244, 722 (2007).
  * Sicardi  _et al._ (2009) E. A. Sicardi, H. Fort, M. H. Vainstein, and J. J. Arenzon, J. Theor. Biol. 256, 240 (2009).
  * Chen  _et al._ (2012) X. Chen, A. Szolnoki, and M. Perc, Phys. Rev. E 86, 036101 (2012).
  * Xiao  _et al._ (2022) S. Xiao, L. Zhang, H. Li, Q. Dai, and J. Yang, Eur. Phys. J. B 95, 67 (2022).
  * Park  _et al._ (2013) J. Park, Y. Do, Z.-G. Huang, and Y.-C. Lai, Chaos 23, 023128 (2013).
  * Wang  _et al._ (2022) X. Wang, Y. Lu, L. Shi, and J. Park, Sci. Rep. 12, 1821 (2022).
  * Park (2019) J. Park, EPL 126, 38004 (2019).
  * Lin  _et al._ (2011) Y.-T. Lin, H.-X. Yang, Z.-X. Wu, and B.-H. Wang, Physica A 390, 77 (2011).
  * Wu  _et al._ (2012) T. Wu, F. Fu, Y. Zhang, and L. Wang, Phys. Rev. E 85, 066104 (2012).
  * Chen  _et al._ (2016a) Y.-S. Chen, H.-X. Yang, and W.-Z. Guo, Physica A 450, 506 (2016a).
  * Cong  _et al._ (2017) R. Cong, Q. Zhao, K. Li, and L. Wang, Sci. Rep. 7, 14015 (2017).
  * Park (2022) J. Park, Chaos, Solit. Fract. 158, 112019 (2022).
  * Li  _et al._ (2020) Q. Li, Z. Zhang, K. Li, L. Chen, Z. Wei, and J. Zhang, Physica A 545, 123664 (2020).
  * Chen  _et al._ (2016b) W. Chen, T. Wu, Z. Li, and L. Wang, Physica A 443, 192 (2016b).
  * Xiao  _et al._ (2020) Z. Xiao, X. Chen, and A. Szolnoki, New J. Phys. 22, 023012 (2020).
  * Dhakal  _et al._ (2022) S. Dhakal, R. Chiong, M. Chica, and T. A. Han, R. Soc. Open Sci. 9, 212000 (2022).
  * Lütz  _et al._ (2021) A. F. Lütz, M. A. Amaral, and L. Wardil, Phys. Rev. E 104, 014304 (2021).
  * Cardinot  _et al._ (2019) M. Cardinot, C. O’Riordan, J. Griffith, and A. Szolnoki, New J. Phys. 21, 073038 (2019).
  * Liu  _et al._ (2012) Y. Liu, X. Chen, L. Zhang, F. Tao, and L. Wang, Chaos Solit. Fract. 45, 1301 (2012).
  * Sigmund (2010) K. Sigmund,  _The Calculus of Selfishness_ (Princeton University Press, Princeton, NJ, 2010).
  * Perc  _et al._ (2013) M. Perc, J. Gómez-Gardeñes, A. Szolnoki, and L. M. Floría and Y. Moreno, J. R. Soc. Interface 10, 20120997 (2013).
  * Wang  _et al._ (2012a) Z. Wang, A. Szolnoki, and M. Perc, Phys. Rev. E 85, 037101 (2012a).
  * Wang  _et al._ (2012b) Z. Wang, A. Szolnoki, and M. Perc, Sci. Rep. 2, 369 (2012b).
  * Szolnoki  _et al._ (2009) A. Szolnoki, M. Perc, and G. Szabó, Phys. Rev. E 80, 056109 (2009).
  * Perc  _et al._ (2008) M. Perc, A. Szolnoki, and G. Szabó, Phys. Rev. E 78, 066101 (2008).
  * Szolnoki and Perc (2009) A. Szolnoki and M. Perc, Eur. Phys. J. B 67, 337 (2009).
  * Wei  _et al._ (2021) X. Wei, P. Xu, S. Du, G. Yan, and H. Pei, Eur. Phys. J. B 94, 210 (2021).
  * Shen  _et al._ (2022) Y. Shen, W. Yin, H. Kang, H. Zhang, and M. Wang, Phys. Lett. A 428, 127935 (2022).
  * Liang  _et al._ (2022) R. Liang, Q. Wang, J. Zhang, G. Zheng, L. Ma, and L. Chen, Phys. Rev. E 105, 054302 (2022).
  * Szolnoki  _et al._ (2011) A. Szolnoki, G. Szabó, and L. Czakó, Phys. Rev. E 84, 046106 (2011).
  * Vicsek and Zafeiris (2012) T. Vicsek and A. Zafeiris, Phys. Rep. 517, 71 (2012).
  * Zheng  _et al._ (2022) J. Zheng, Y. He, T. Ren, and Y. Huang, Physica A 596, 127101 (2022).
  * Flores  _et al._ (2021) L. S. Flores, H. C. Fernandes, M. A. Amaral, and M. H. Vainstein, J. Theor. Biol. 524, 110737 (2021).
  * Quan  _et al._ (2021a) J. Quan, C. Tang, and X. Wang, Physica A 563, 125488 (2021a).
  * Deng and Zhang (2022) Y. Deng and J. Zhang, Eur. Phys. J. B 95, 29 (2022).
  * Lee  _et al._ (2022) H.-W. Lee, C. Cleveland, and A. Szolnoki, Applied Mathematics and Computation 417, 126797 (2022).
  * Liu  _et al._ (2021) J. Liu, M. Peng, Y. Peng, Y. Li, C. Chu, X. Li, and Q. Liu, Eur. Phys. J. B 94, 167 (2021).
  * Lv and Song (2022) S. Lv and F. Song, Appl. Math. Comput. 412, 126586 (2022).
  * Quan  _et al._ (2021b) J. Quan, Z. Pu, and X. Wang, Chaos, Solitons and Fractals 151, 111229 (2021b).
  * Zheng  _et al._ (2021) L. Zheng, H. Xu, C. Tian, and S. Fan, Physica A 581, 126228 (2021).
  * Li  _et al._ (2021a) X. Li, L. Cheng, X. Niu, S. Li, C. Liu, and P. Zhu, Eur. Phys. J. B 94, 126 (2021a).
  * Fu  _et al._ (2021) M. Fu, J. Wang, L. Cheng, and L. Chen, Physica A 580, 125672 (2021).
  * Zhu  _et al._ (2021) P. Zhu, X. Hou, Y. Guo, J. Xu, and J. Liu, Eur. Phys. J. B 94, 58 (2021).
  * Yang and Zhang (2021) L. Yang and L. Zhang, Chaos Solit. and Fract. 142, 110485 (2021).
  * Li  _et al._ (2021b) K. Li, Y. Mao, Z. Wei, and R. Cong, Chaos, Solit. and Fract. 143, 110591 (2021b).
  * Lee  _et al._ (2021) H.-W. Lee, C. Cleveland, and A. Szolnoki, Physica A 582, 126222 (2021).
  * Quan  _et al._ (2020) J. Quan, Y. Qin, Y. Zhou, X. Wang, and J.-B. Yang, J. Stat. Mech. 2020, 093405 (2020).
  * Fu  _et al._ (2019) M. Fu, W. Guo, L. Cheng, S. Huang, and D. Chen, Physica A 525, 1323 (2019).