# Nonlinear molecular excitations in a completely inhomogeneous DNA chain

M. Daniel<EMAIL_ADDRESS>V. Vasumathi Corresponding Author.
Telephone:+91-431-2407057, Fax:+91-431-2407093  Centre for Nonlinear Dynamics,
School of Physics, Bharathidasan University, Tiruchirappalli - 620 024, India.

###### Abstract

We study the nonlinear dynamics of a completely inhomogeneous DNA chain which
is governed by a perturbed sine-Gordon equation. A multiple scale perturbation
analysis provides perturbed kink-antikink solitons to represent open state
configuration with small fluctuation. The perturbation due to inhomogeneities
changes the velocity of the soliton. However, the width of the soliton remains
constant.

###### keywords:

DNA , Soliton , Multiple Scale Perturbation

###### PACS:

87.15.He, 66.90.+r, 63.20.Ry

,

##  1 Introduction

A number of theoretical models have been proposed in the recent times to study
the nonlinear dynamics of Deoxyribonucleic acid (DNA) molecule to understand
the conservation and transformation of genetic information (see for e.g [1,
2]). These models are based on longitudinal, transverse, rotational and
stretching motions of bases. Among these different possible motions,
rotational motion of bases is found to contribute more towards the opening of
base pairs and to the nonlinear dynamics of DNA. The first contribution
towards nonlinear dynamics of DNA was made by Englander and his co-workers [3]
and they studied the base pair opening in DNA by taking into account the
rotational motion. Yomosa [4, 5] proposed a plane base rotator model by taking
into account the rotational motion of bases in a plane normal to the helical
axis, and Takeno and Homma generalized the same [6, 7, 8]. Later using this
model, several authors found solitons to govern the fluctuation of DNA double
helix between an open state and its equilibrium states [9, 10, 11, 12, 13,
14]. Peyrard and Bishop [15, 16] and Christiansen and his collegues [17]
studied the process of base pair opening by taking into account the transverse
and longitudinal motions of bases in DNA. Very recently, there have been
extensions of the radial model of Bishop and Peyrard [18, 19], composite
models for DNA torsion dynamics [20] and models interplaying between radial
and torsional dynamics [21, 22, 23, 24]. In all the above studies, homogeneous
strands and hydrogen bonds have been considered for the analysis.

However, in nature the presence of different sites along the strands such as
promotor, coding, terminator, etc., each of which has a specific sequence of
bases is related to a particular function and thus making the strands site-
dependent or inhomogeneous [25, 26]. Also, the presence of abasic sites leads
to inhomogeneity in stacking [27]. In this context, in a recent paper the
present authors [28] studied the nonlinear molecular excitations in DNA with
site-dependent stacking energy along the strands based on the plane base
rotator model. The nonlinear dynamics of DNA in this case was found to be
governed by a perturbed sine-Gordon (s-G) equation. The perturbed kink and
antikink soliton solutions of the perturbed s-G equation represented an open
state configuration of base pairs with small fluctuation. The perturbation in
this case introduces small fluctuations in the localized region of the soliton
retaining the overall shape of the soliton. However, the width of the soliton
remains constant and the velocity changes for different inhomogeneities. The
results indicate that the presence of inhomogeneity in stacking changes the
number of base pairs that participate in the open state configuration and
modifies the speed with which the open state configuration travels along the
double helical chain. In reality, the presence of site-dependent strands in
DNA changes the nature of hydrogen bonds between adjacent base pairs and the
presence of abasic sites leads to absence of hydrogen bonds. Thus, when the
strands are site-dependent in stacking, naturally the hydrogen bonds that
connect the bases between the strands are also site-dependent. Hence, it has
become necessary to consider inhomogeneity in hydrogen bonds also in the study
of nonlinear dynamics of DNA. In the present paper, we study the dynamics of
DNA with inhomogeneity both in stacking and in hydrogen bonds using the plane
base rotator model. The paper is organized as follows. In section 2, we
present the Hamiltonian of our model and derive the associated dynamical
equation for the inhomogeneous DNA. The effect of stacking and hydrogen bond
inhomogeneity on base-pair opening is studied by solving the dynamical
equations using a multiple scale soliton perturbation theory in section 3. The
results are concluded in section 4.

##  2  Hamiltonian and the dynamical equation

We consider the B-form of a DNA double helix with site-dependent strands as
well as base-pair sequence and study the nonlinear molecular excitations by
considering a plane-base rotator model. In Fig. (1a) we have presented a
sketch of the DNA double helix. Here, $S$ and $S^{\prime}$ represent the two
complementary strands in the DNA double helix and each arrow represents the
direction of the base attached to the strand and the dots between arrows
represent the net hydrogen bonding effect between the complementary bases. The
z-axis is chosen along the helical axis of the DNA. Fig. (1b) represents a
horizontal projection of the $n^{th}$ base pair in the xy-plane. In this
figure $Q_{n}$ and $Q^{\prime}_{n}$ denote the tips of the $n^{th}$ bases
belonging to the strands $S$ and $S^{\prime}$. $P_{n}$ and $P^{\prime}_{n}$
represent the points where the bases in the $n^{th}$ base pair are attached to
the strands $S$ and $S^{\prime}$ respectively.

As we are looking for opening of base pairs in DNA which is related to
important DNA functions such as replication and transcription, we consider the
rotational motion of bases (due to its importance among other motions) in a
plane normal to the helical axis (z-direction) represented by the angles
$\phi_{n}$ and $\phi^{\prime}_{n}$ at the $n^{th}$ site of the base pair. The
stacking and hydrogen bonding energies are the major components of the energy
in a DNA double helix. In the case of a homogeneous DNA system, Yomosa [4, 5]
expressed the Hamiltonian involving these energies in terms of the rotational
angles $\phi_{n}$ and $\phi^{\prime}_{n}$ under the plane base rotator model
which was later modified by the present authors [28] in the case of site-
dependent stacking. When both the stacking and hydrogen bonds are site-
dependent, the Hamiltonian for our plane base rotator model of DNA double
helix is written in terms of the rotational angles as

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{n}\left[\frac{I}{2}({\dot{\phi}_{n}}^{2}+{\dot{\phi}_{n}}^{{}^{\prime}2})+Jf_{n}\left[2-\cos(\phi_{n+1}-\phi_{n})-\cos(\phi^{\prime}_{n+1}-\phi^{\prime}_{n})\right]\right.$
(1) $\displaystyle\left.-\eta
g_{n}\left[1-\cos(\phi_{n}-\phi^{\prime}_{n})\right]\right].$

The first two terms in the Hamiltonian (1) represent the kinetic energies of
the rotational motion of the $n^{th}$ nucleotide bases with $I$ their moments
of inertia and the remaining terms represent the potential energy due to
stacking and hydrogen bonds. While $J$ and $\eta$ represent a measure of
stacking and hydrogen bonding energies respectively, $f_{n}$ and $g_{n}$
indicate the site-dependent(inhomogeneous) character of stacking and hydrogen
bonds respectively.
The Hamilton’s equations of motion corresponding to Hamiltonian (1) is written
as

$\displaystyle I\frac{\partial^{2}\phi_{n}}{\partial t^{2}}$ $\displaystyle=$
$\displaystyle
J\left[f_{n}\sin(\phi_{n+1}-\phi_{n})-f_{n-1}\sin(\phi_{n}-\phi_{n-1})\right]$
(2a) $\displaystyle+\eta g_{n}\sin(\phi_{n}-\phi^{\prime}_{n}),$
$\displaystyle I\frac{\partial^{2}\phi^{\prime}_{n}}{\partial t^{2}}$
$\displaystyle=$ $\displaystyle
J\left[f_{n}\sin(\phi^{\prime}_{n+1}-\phi^{\prime}_{n})-f_{n-1}\sin(\phi^{\prime}_{n}-\phi^{\prime}_{n-1})\right]$
(2b) $\displaystyle+\eta g_{n}\sin(\phi^{\prime}_{n}-\phi_{n}).$

It is expected that the difference in the angular rotation of neighbouring
bases along the two strands in the case of B-form of DNA double helix is small
[6, 7]. Therefore under small angle approximation, in the continuum limit Eqs.
(2a) and (2b) are written as

$\displaystyle\phi_{\hat{t}\hat{t}}$ $\displaystyle=$ $\displaystyle
f(z)\phi_{zz}+f_{z}\phi_{z}-\frac{1}{2}g(z)\sin(\phi-\phi^{\prime}),$ (3a)
$\displaystyle\phi^{\prime}_{\hat{t}\hat{t}}$ $\displaystyle=$ $\displaystyle
f(z)\phi^{\prime}_{zz}+f_{z}\phi^{\prime}_{z}-\frac{1}{2}g(z)\sin(\phi^{\prime}-\phi).$
(3b)

While writing Eqs. (3a) and (3b) we have chosen $\eta=-\frac{Ja^{2}}{2}$ and
rescaled the time variable as $\hat{t}=\sqrt{\frac{Ja^{2}}{I}}t$. Now, adding
and subtracting Eqs.(3a) and (3b) and by choosing $\phi=-\phi^{\prime}$, we
obtain

$\displaystyle\Psi_{\hat{t}\hat{t}}-\Psi_{zz}+\sin\Psi=\epsilon\left[A\left(s(z)\Psi_{z}\right)_{z}-B~{}h(z)\sin\Psi\right],$
(4)

where $\Psi=2\phi$. Further, as the contribution due to inhomogeneity is small
compared to homogeneous stacking and hydrogen bonding energies while writing
Eq.(4), the inhomogeneity in stacking and hydrogen bonding energies are
expressed in terms of a small parameter $\epsilon$ as
$f(z)=1+\epsilon~{}As(z)$ and $g(z)=1+\epsilon~{}Bh(z)$, where $A$ and $B$ are
arbitrary constants. When $\epsilon=0,$ Eq. (4) reduces to the completely
integrable sine-Gordon(s-G) equation which admits N-soliton solutions in the
form of kink and antikink [29]. Hence, we call Eq. (4) as a perturbed sine-
Gordon equation. For instance, the one soliton solution of the integrable s-G
equation obtained through Inverse Scattering Transform (IST) method is written
as,

$\displaystyle\Psi(z,\hat{t})=4arctan{exp[\pm m_{0}(z-v_{0}\hat{t})]},\quad
m_{0}=\frac{1}{\sqrt{1-v_{0}^{2}}}.$ (5)

Here $v_{0}$ and $m^{-1}_{0}$ are real parameters that represent the velocity
and width of the soliton respectively. In Eq.(5) while the upper sign
corresponds to the kink soliton, the lower sign represents the antikink
soliton. In Figs. (2a) and (2b) the one soliton kink-antikink solutions
(Eq.(5)) are plotted by choosing $v_{0}=0.4$. The kink-antikink soliton
solution of the integrable sine-Gordon equation describes an open state
configuration in the DNA double helix. In Figure (2c), we present a sketch of
how the base pairs open locally in the form of kink-antikink structure in each
strand and propagate along the direction of the helical axis.

The base pair opening will help in the process of replication that duplicates
DNA and in transcription which helps to synthesize messenger RNA. However,
when $\epsilon\neq 0$, the inhomogeneity in stacking and hydrogen bonds may
affect the base-pair opening through a perturbation on the kink-antikink
solitons of the s-G equation. Therefore in the next section, we solve the
perturbed s-G equation (4) using a multiple-scale soliton perturbation theory
[30, 31, 32] (as has been carried out in the case corresponding to $B=0$ [28])
to understand the effect of stacking and hydrogen bond inhomogeneities on the
base pair opening.

##  3 Effect of stacking and hydrogen bonding inhomogeneities on the open
state configuration

###  3.1 Multiple-scale soliton perturbation theory

In order to study the effect of inhomogeneity in stacking and hydrogen bonds
on the base pair opening in the form of kink-antikink soliton by treating them
as a perturbation, the time variable $\hat{t}$ is transformed into several
variables as $t_{n}=\epsilon^{n}\hat{t}$ where $n=0,1,2,...$ and $\epsilon$ is
a very small parameter [30, 31, 32]. In view of this, the time derivative and
$\Psi$ in Eq. (4) are replaced by the expansions
$\frac{\partial}{\partial\hat{t}}=\frac{\partial}{\partial
t_{0}}+\epsilon~{}\frac{\partial}{\partial
t_{1}}+\epsilon^{2}\frac{\partial}{\partial t_{2}}+...$ and
$\Psi=\Psi^{(0)}+\epsilon\Psi^{(1)}+\epsilon^{2}\Psi^{(2)}+....$ We then
equate the coefficients of different powers of $\epsilon$ and obtain the
following equations. At $O(\epsilon^{(0)})$ we obtain the integrable sine-
Gordon equation

$\displaystyle\Psi^{(0)}_{t_{0}t_{0}}-\Psi^{(0)}_{zz}+sin\Psi^{(0)}=0,$ (6)

for which the one soliton solution takes the form as given in Eq.(5) with
$\hat{t}$ replaced by $t_{0}$. Due to perturbation, the soliton parameters
namely $m$ and $\xi(\xi=v\hat{t})$ are now treated as functions of the slow
time variables $t_{0},t_{1},t_{2},....$ However, $m$ is treated as independent
of $t_{0}$. The equation at $O(\epsilon^{(1)})$ takes the form

$\displaystyle{\Psi^{(1)}_{\tau\zeta}}-{\Psi_{\zeta\zeta}^{(1)}}+(1-2{sech^{2}}\zeta)\Psi^{(1)}=F^{(1)}(\zeta,\tau),$
(7) $\displaystyle
where\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad$
$\displaystyle
F^{(1)}(\zeta,t_{0})=2\left[As(\zeta)sech\zeta\right]_{\zeta}-2aBh(\zeta)\tanh\zeta\mbox{sech}\zeta\qquad\qquad$
$\displaystyle+4v_{0}sech\zeta\left[m_{t_{1}}+(m^{2}\xi_{t_{1}}-\zeta
m_{t_{1}})\tanh\zeta\right].$ (8)

While writing the above equation we have used the transformation
$\hat{\zeta}=m(z-vt_{0})$ and $\hat{t}_{0}=t_{0}$ to represent everything in a
co-ordinate frame moving with the soliton. Then, we have used another set of
transformations given by
$\tau=\frac{\hat{t}_{0}}{2m}-\frac{(1+v)\hat{\zeta}}{2}$ and
$\zeta=\hat{\zeta}$ for our later convenience. We have also replaced
$sin\Psi^{(0)}$ by $2a\tanh\zeta$ $sech\zeta$, where $a=\pm 1$. The solution
of Eq. (7) is searched by assuming $\Psi^{(1)}(\zeta,\tau)=X(\zeta)T(\tau)$
and $F^{(1)}(\zeta,\tau)=X_{\zeta}(\zeta)H(\tau).$ Substituting the above in
Eq. (7) and simplifying we obtain

$\displaystyle X_{\zeta\zeta}+(2sech^{2}\zeta-1)X=\lambda_{0}X_{\zeta},$ (9)
$\displaystyle T_{\tau}-\lambda_{0}T=H(\tau),$ (10)

where $\lambda_{0}$ is a constant. Thus, the problem of constructing the
perturbed soliton at this moment turns out to be solving Eqs. (9) and (10) by
constructing the eigenfunctions and finding the eigen values [28, 32].
It may be noted that Eq.(9) differs from the normal eigen value problem, with
$X_{\zeta}$ in the right hand side instead of $X$. Hence, before actually
solving the eigen value equation (9), we first consider it in a more general
form given by

$\displaystyle L_{1}X=\lambda\tilde{X},\qquad
L_{1}=\partial_{\zeta\zeta}+2sech^{2}\zeta-1,$ (11)

where $\lambda$ is the eigen value. To find the adjoint eigen function to $X$,
we consider another eigen value problem

$\displaystyle L_{2}\tilde{X}=\lambda X,$ (12)

where $L_{2}$ is to be determined. Combining the two eigen value problems we
get

$\displaystyle
L_{2}L_{1}X=\lambda^{2}X,~{}~{}L_{1}L_{2}\tilde{X}=\lambda^{2}\tilde{X}.$ (13)

From the above equations we conclude that the operator $L_{1}L_{2}$ is the
adjoint of $L_{2}L_{1}$ and also $X$ and $\tilde{X}$ are expected to be
adjoint eigen functions. Hence, we can find the eigenfunction by solving
Eqs.(11) and (12). However, eventhough $L_{1}$ is known as given in Eq.(11),
the operator $L_{2}$ is still unknown. So, by experience we choose
$L_{2}=\partial_{\zeta\zeta}+6sech^{2}\zeta-1$, and solve the eigen value
equations by choosing the eigen functions as

$\displaystyle
X(\zeta,k)=p(\zeta,k)e^{ik\zeta},~{}~{}\tilde{X}(\zeta,k)=q(\zeta,k)e^{ik\zeta},$
(14)

where $k$ is the propagation constant. Further, as the operator $L_{2}L_{1}$
is self-adjoint, non-negative and satisfying a regular eigenvalue problem, the
sine-Gordon soliton is expected to be stable. On substituting Eqs. (14) in
Eqs.(11) and (12) in the asymptotic limit, we obtain the eigen value as
$\lambda=-(1+k^{2})$. Now, in order to find the eigen functions we expand
$p(\zeta,k)$ and $q(\zeta,k)$ [30, 32] as

$\displaystyle p(\zeta,k)$ $\displaystyle=$ $\displaystyle
p_{0}+p_{1}\frac{\sinh\zeta}{\cosh\zeta}+p_{2}\frac{1}{\cosh^{2}\zeta}+p_{3}\frac{\sinh\zeta}{\cosh^{3}\zeta}+p_{4}\frac{1}{\cosh^{4}\zeta}+...,$
(15a) $\displaystyle q(\zeta,k)$ $\displaystyle=$ $\displaystyle
q_{0}+q_{1}\frac{\sinh\zeta}{\cosh\zeta}+q_{2}\frac{1}{\cosh^{2}\zeta}+q_{3}\frac{\sinh\zeta}{\cosh^{3}\zeta}+q_{4}\frac{1}{\cosh^{4}\zeta}+...,$
(15b)

where $p_{j}$ and $q_{j}$, j=0,1,2,… are functions of $k$ to be determined. On
substituting Eqs. (14), (15a) and (15b) in Eqs. (11) and (12) and collecting
the coefficients of $1,~{}\frac{\sinh\zeta}{\cosh\zeta}$,
$\frac{1}{\cosh^{2}\zeta}$,… we get a set of simultaneous equations. On
solving those equations we obtain the following eigen functions [28, 32].

$\displaystyle X(\zeta,k)$ $\displaystyle=$
$\displaystyle\frac{(1-k^{2}-2ik\tanh\zeta)}{\sqrt{2\pi}(1+k^{2})}e^{ik\zeta},$
(16) $\displaystyle\tilde{X}(\zeta,k)$ $\displaystyle=$
$\displaystyle\frac{(1-k^{2}-2ik\tanh\zeta-2sech^{2}\zeta)}{\sqrt{2\pi}(1+k^{2})}e^{ik\zeta}.$
(17)

The higher order coefficients $p_{3},p_{4},...$ and $q_{3},q_{4},...$ vanish.
It may be noted that Eq.(10) is a linear inhomogeneous differential equation
and can be solved using known procedures. The solution reads

$\displaystyle
T(\tau,k)=\frac{1}{i\lambda_{0}k(1+k^{2})}\int_{-\infty}^{\infty}d\zeta^{\prime}F^{(1)}(\zeta^{\prime},\tau){X}^{\ast}(\zeta^{\prime},k)(e^{\lambda_{0}[\tau-\frac{(1+v)}{2}\zeta^{\prime}]}-1),$
(18)

where $\lambda_{0}=\frac{i(1+k^{2})}{k}$. The first order correction to the
soliton can be computed using the following expression.

$\displaystyle\Psi^{(1)}(\zeta,\tau)=\int_{-\infty}^{\infty}X(\zeta,k)T(\tau,k)dk+\sum_{j=0,1}X_{j}(\zeta)T_{j}(\tau).$
(19)

Here the continuous eigenfunctions $X(\zeta,k)$ and $T(\tau,k)$ are already
known as given in Eqs. (16) and (18). However, the discrete eigen states
$X_{0},X_{1}$ and $T_{0},T_{1}$ are unknown. The two discrete eigenstates
$X_{0},X_{1}$ corresponding to the discrete eigen value $\lambda=0$ can be
found using the completeness of the continuous eigenfunctions as

$\displaystyle X_{0}(\zeta)$ $\displaystyle=$ $\displaystyle
sech\zeta,~{}X_{1}(\zeta)=\zeta sech\zeta.$ (20)

In order to find $T_{0}$ and $T_{1}$, we substitute Eq.(19) in Eq.(7) and
multiply by $X_{0}(\zeta)$ and $X_{1}(\zeta)$ separately and use the
orthonormal relations to get

$\displaystyle{T_{1}}_{\tau}(\tau)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}F^{(1)}(\zeta,\tau)X_{0}(\zeta)d\zeta,$
(21a) $\displaystyle{T_{0}}_{\tau}(\tau)-2T_{1}(\tau)$ $\displaystyle=$
$\displaystyle-\int_{-\infty}^{\infty}F^{(1)}(\zeta,\tau)X_{1}(\zeta)d\zeta.\qquad$
(21b)

As $F^{(1)}(\zeta,\tau)$ given in Eq.(8) does not contain the time variable
$\tau$ explicitly, the right hand side of Eqs. (21a) and (21b) also should be
independent of time, and hence we write the nonsecularity conditions as

$\displaystyle\int_{-\infty}^{\infty}F^{(1)}(\zeta,\tau)X_{0}(\zeta)d\zeta=0,$
(22a)
$\displaystyle\int_{-\infty}^{\infty}F^{(1)}(\zeta,\tau)X_{1}(\zeta)d\zeta=0.$
(22b)

It may be indicated that, the nonsecularity conditions give way for the
stability of soliton [33] under perturbation through an understanding of the
evolution of soliton parameters such as width and velocity. Substituting the
above equations in Eqs. (21a) and (21b), we choose $T_{1}(\tau)=0$ and find
$T_{0}(\tau)=C$ which has to be determined. For this, we substitute
$T_{1}(\tau)=0$ in Eq. (21b) and integrate to obtain

$\displaystyle T_{0}(\tau)$ $\displaystyle=$
$\displaystyle\frac{(1+v)}{2}\int_{-\infty}^{\infty}d\zeta~{}\zeta
F^{(1)}(\zeta,\tau)X_{1}(\zeta).$ (23)

In order to find the first order correction to soliton $\Psi^{(1)}$ we need to
evaluate the eigen states $T(\tau,k)$ and $T_{0}(\tau)$ explicitly for which
we have to find the explicit form of $F^{(1)}(\zeta,\tau)$ which contains
unknown quantities like $m_{t_{1}},~{}\xi_{t_{1}},~{}s(\zeta)$ and $h(\zeta)$.
Therefore, we substitute Eqs. (8) and (20) in the nonsecularity conditions
given in Eqs. (22a) and (22b) and evaluate the integrals to obtain

$\displaystyle m_{t_{1}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2v_{0}}\int_{-\infty}^{\infty}\left(A\left[{s(\zeta)sech\zeta}\right]_{\zeta}-aBh(\zeta)sech\zeta\tanh\zeta\right)sech\zeta
d\zeta,$ (24a) $\displaystyle\xi_{t_{1}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2m^{2}v_{0}}\int_{-\infty}^{\infty}\left(A\left[{s(\zeta)sech\zeta}\right]_{\zeta}-aBh(\zeta)sech\zeta\tanh\zeta\right)\zeta
sech\zeta d\zeta.$ (24b)

Here $m_{t_{1}}$ and $\xi_{t_{1}}$ represent the time variation of the inverse
of the width and velocity of the soliton respectively.

###  3.2  Variation of soliton parameters

In order to find the variation of the width and velocity of the soliton while
propagating along the inhomogeneous DNA chain, we assume that the soliton with
a width of $m_{0}^{-1}$ is travelling with a speed of $v_{0}(\xi_{t_{0}})$
when the perturbation is switched on. In other words, the problem now boils
down to understand the propagation of the kink-antikink soliton in an
inhomogeneous DNA chain and to measure the change of the soliton parameters
due to inhomogeneity. To find the variation of the soliton parameters
explicitly and to construct the perturbed soliton solution we have to evaluate
the integrals found in the right hand sides of Eqs. (24a) and (24b) which can
be carried out only on supplying specific forms of $s(\zeta)$ and $h(\zeta)$.
Hence, we consider inhomogeneity in the form of localized and periodic
functions separately and further assume that the inhomogeneity in the stacking
and in the hydrogen bonds are equal. The localized inhomogeneity represents
the intercalation of a compound between neighbouring base pairs or the
presence of defect or the presence of abasic site in the DNA double helical
chain. The periodic nature of inhomogeneity may represent a periodic
repetition of similar base pairs along the helical chain.

####  3.2.1  Localized inhomogeneity

To understand the effect of localized inhomogeneity in stacking and hydrogen
bonds on the width and velocity of the soliton during propagation, we
substitute $s(\zeta)=h(\zeta)=sech\zeta$ in Eqs.(24a) and (24b). On evaluating
the integrals, we obtain
$m_{t_{1}}=0,~{}\xi_{t_{1}}=\frac{\pi}{12m^{2}v_{0}}(2A+aB)$ which can be
written in terms of the original time variable $\hat{t}$ by using the
expansions $m_{\hat{t}}=m_{t_{0}}+\epsilon m_{t_{1}}$ and
$\xi_{\hat{t}}=\xi_{t_{0}}+\epsilon\xi_{t_{1}}$ as

$\displaystyle m=m_{0},\quad\xi_{\hat{t}}\equiv
v=v_{0}+\frac{\epsilon\pi(2A+aB)}{12m^{2}v_{0}},$ (25)

where $1/m_{0}$ is the width and $v_{0}$ is the velocity of the soliton in the
absence of inhomogeneity. The first of Eq.(25) says that when the
inhomogeneities in stacking and in hydrogen bonds are in the form
$s(\zeta)=h(\zeta)=sech\zeta$, the width of the soliton remains constant,
thereby showing that the number of base pairs participating in the opening
process remain constant during propagation. However, from the second of
Eq.(25), we find that the velocity of the soliton gets a correction. The
nature of correction in the velocity depends on the value of ‘$a$’ which takes
$+1$ or $-1$ and also on the nature of $A$ and $B$ which can be either
positive or negative. When $A$ and $B$ are greater than zero, the
inhomogeneity will correspond to an energetic barrier and on the other hand
when $A$ and $B$ are less than zero, the inhomogeneity will correspond to $a$
potential well. First, we consider the case corresponding to $a=+1$. In this
case when $(2A+B)>0$, the velocity of the soliton gets a positive correction
and it may propagate along the chain without formation of a bound state. On
the other hand, when $(2A+B)<0$, the inhomogeneities slow down the soliton.
Ofcourse, when $(2A+B)=0$, that is when the inhomogeneities in the stacking
and in the hydrogen bonding suitably balance each other, the velocity of the
soliton remains unaltered. Finally, the soliton stops when the original
velocity satisfies the condition
$v_{0}^{2}=\frac{1}{1-\frac{12}{\epsilon\pi(2A+B)}}$. In all the above cases,
the stability of the soliton is guaranteed. A similar argument can be made in
the case when $a=-1$ with $(2A+B)$ replaced by $(2A-B)$. It may also be noted
that similar results have been obtained in the case of resonant kink-impurity
interaction and kink scattering in a perturbed sine-Gordon model by Zhang Fei
et al [34], say that the kink will pass the impurity and escape to the
positive infinity when the initial velocity of kink soliton is larger than the
critical value. In a similar context Yakushevich et al [35] while studying the
interaction between soliton and the point defect in DNA chain showed
numerically, that the solitons are stable. At this point it is also worth
mentioning that Dandoloff and Saxena [36] realized that in the case of an XY-
coupled spin chain model which is identifiable with our plane-base rotator
model of DNA, the ansatz $sech\zeta$ energetically favours the deformation of
spin chain.

####  3.2.2  Periodic inhomogeneity

We then choose the periodic inhomogeneity in the form
$s(\zeta)=h(\zeta)=\cos\zeta$ and substitute the same in Eqs. (24a) and (24b).
On evaluating the integrals, we get $m_{t_{1}}=0,~{}\xi_{t_{1}}=\frac{\pi[\pi
A+a(4-\pi)B]}{16m^{2}v_{0}}$ which can be written in terms of the original
time variable $\hat{t}$ as

$\displaystyle m=m_{0},\quad\xi_{\hat{t}}\equiv v=v_{0}+\frac{\epsilon\pi[\pi
A+a(4-\pi)B]}{16m^{2}v_{0}}.$ (26)

From Eq.(26), we observe that the width of the soliton remains constant and
the velocity gets a correction. Here also, one can project an argument similar
to the case of localized inhomogeneity. The only difference between the two
cases is the quantum of correction added to soliton velocity. Normally in
DNAs, inhomogeneity in hydrogen bonds is expected to be dominant and therefore
one would expect that $B>A$. One can verify that it is possible to obtain the
above condition from the velocity corrections in Eqs. (25) and (26) by writing
$\frac{\epsilon\pi[\pi
A+(4-\pi)B]}{16m^{2}v_{0}}>\frac{\epsilon\pi(2A+B)}{12m^{2}v_{0}}$ when
$a=+1$. The above condition indicates that the correction in velocity in the
case of periodic inhomogeneity is larger than that in the case of localized
inhomogeneity. This is because in this case, the inhomogeneity occurs
periodically in the entire length of the DNA chain. In a recent paper,
Yakushevich et al [35] studied numerically the dynamics of topological
solitons describing open states in an inhomogeneous DNA and investigated
interaction of soliton with the inhomogeneity and the results, have very close
analogy with the results of our perturbation analysis. It was shown that the
soliton can easily propagate along the DNA chain without forming a bound state
thus showing that soliton moving with sufficiently large velocity along the
DNA chain is stable with respect to defect or inhomogeneity.

###  3.3 First order perturbed soliton

Having understood the variation of the width and velocity of the soliton in a
slow time scale due to perturbation, we now construct the first order
perturbed soliton by substituting the values of the basis functions
$\\{X\\}\equiv\\{X(\zeta,k),X_{0}(\zeta),X_{1}(\zeta)\\}$ and
$\\{T\\}\equiv\\{T(\tau,k),T_{0}(\tau),T_{1}(\tau)\\}$ given in Eqs. (16),
(20), (18), (23) with $T_{1}(\tau)=0$ in Eq.(19), we get

$\displaystyle\Psi^{(1)}(\zeta,t_{0})$ $\displaystyle=$
$\displaystyle-\frac{1}{\pi}\int_{-\infty}^{\infty}dk\frac{1}{(1+k^{2})^{3}}(1-k^{2}-2ik\tanh\zeta)e^{ik\zeta}$
(27)
$\displaystyle\times\int_{-\infty}^{\infty}d\zeta^{\prime}\left(\left[As(\zeta^{\prime})sech\zeta^{\prime}\right]_{\zeta^{\prime}}-aBh(\zeta^{\prime})sech\zeta^{\prime}\tanh\zeta^{\prime}\right.+2v_{0}~{}sech\zeta^{\prime}$
$\displaystyle\times\left.\left[m_{t_{1}}+(m^{2}\xi_{t_{1}}-\zeta^{\prime}m_{t_{1}})\tanh\zeta^{\prime}\right]\right)(1-k^{2}+2ik\tanh\zeta^{\prime})e^{-ik\zeta^{\prime}}$
$\displaystyle\times\left(e^{\lambda_{0}[\tau+\frac{(1+v_{0})}{2}\zeta^{\prime}]}-1\right)+(1+v)sech\zeta\int_{-\infty}^{\infty}d\zeta^{\prime}\left(\left[As(\zeta^{\prime})sech\zeta^{\prime}\right]_{\zeta^{\prime}}\right.$
$\displaystyle-
aBh(\zeta^{\prime})sech\zeta^{\prime}\tanh\zeta^{\prime}+2v_{0}~{}sech\zeta^{\prime}\left[m_{t_{1}}\right.$
$\displaystyle\left.+(m^{2}\xi_{t_{1}}-\zeta^{\prime}\left.m_{t_{1}})\tanh\zeta^{\prime}\right]\right)\zeta^{\prime
2}sech\zeta^{\prime}.$

While writing the above, we have also used
$\tau=\frac{1}{2m}[t_{0}-m(1+v)\zeta]$. It may be noted that majority of the
poles that lie within the contour in Eq. (27), are purely imaginary giving
rise to exponentially localized residues and hence do not give rise to any
radiation thus will lead to stable form of soliton [37].

####  3.3.1 Localized inhomogeneity

Now, we explicitly construct the first order perturbation correction to the
one soliton in the case of the localized inhomogeneity (i.e) when
$s(\zeta)=h(\zeta)=sech\zeta$, by substituting the corresponding values of
$F^{(1)},m_{t_{1}}$ and $\xi_{t_{1}}$ in Eq. (27). The resultant integrals are
then evaluated using standard residue theorem [38] which involves very lengthy
algebra and at the end we obtain

$\displaystyle\Psi^{(1)}(\zeta,t_{0})$ $\displaystyle\approx$
$\displaystyle(2A+aB)\left[\frac{40}{27\sqrt{2}}\sqrt{sech\zeta}~{}e^{-\frac{3\zeta}{2}}+\frac{32}{27\sqrt[3]{2}}\tanh\zeta\sqrt[3]{sech\zeta}~{}e^{-\frac{5}{3}\zeta}\right.$
(28)
$\displaystyle+\frac{\pi}{12v^{2}}\left.[2vt_{0}+m^{2}(v^{2}-1)]sech\zeta\right].$

Finally, the perturbed one soliton solution, that is
$\Psi(z,t_{0})=\Psi^{(0)}(z,t_{0})+\Psi^{(1)}(z,t_{0})$ (choosing
$\epsilon=1$) is written in terms of the original variables as

$\displaystyle\Psi(z,t_{0})$ $\displaystyle\approx$ $\displaystyle
4arc\tan\exp[\pm
m_{0}(z-v_{0}t_{0})]+(2A+aB)\left[\frac{40}{27\sqrt{2}}e^{\mp\frac{3(m(z-vt_{0})}{2}}\right.$
(29) $\displaystyle\times\sqrt{sech[\pm
m(z-vt_{0})]}+\frac{32}{27\sqrt[3]{2}}\tanh[\pm m(z-vt_{0})]$
$\displaystyle\times\sqrt[3]{sech[\pm
m(z-vt_{0})]}e^{\mp\frac{5}{3}m(z-vt_{0})}+\frac{\pi}{12mv^{2}}$
$\displaystyle\left.\times\left[m(v^{2}-1)+2t_{0}v\right]sech[\pm
m(z-vt_{0})]\right].$

Knowing $\Psi$, the angle of rotation of bases $\phi(z,t_{0})$ can be
immediately written down by using the relation $\phi=\frac{\Psi}{2}$. In Figs.
3(a,b) we have plotted $\phi(z,t_{0})$, the rotation of bases under the
perturbation $g(z)=h(z)=sechz$ by choosing $a=A=B=1$ and $v_{0}=0.4$. From the
figures we observe that there appears fluctuation in the form of a train of
pulses closely resembling the shape of the inhomogeneity profile in the width
of the soliton as time progresses. But, there is no change in the topological
character and no fluctuations appear in the asymptotic region of the soliton.
In Figs. (4a) and (4b), we have plotted the perturbed kink and antikink
solitons respectively corresponding to the case when the hydrogen bond
inhomogeneity is absent (B=0) for comparison. From Figs. (3) and (4), we
observe that the fluctuation in the form of a train of pulses appear in both
the cases. Eventhough the perturbed solitons in both the cases appear
qualitatively the same, the inhomogeneity in hydrogen bonds adds more
fluctuation in the width of the soliton.

####  3.3.2 Periodic inhomogeneity

We then repeat the procedure for constructing the perturbed one soliton
solution $\Psi(\zeta,t_{0})$ in the case of periodic inhomogeneity by choosing
$s(\zeta)=h(\zeta)=\cos\zeta$, and obtain the perturbed soliton solution as

$\displaystyle\Psi(z,t_{0})$ $\displaystyle\approx$ $\displaystyle
4~{}arctan\exp[\pm m_{0}(z-v_{0}t_{0})]+[A\pi+aB(4-\pi)][b+dt_{0}]$ (30)
$\displaystyle sech[m_{0}(z-vt_{0})],$

where $b=\frac{\pi^{2}}{16v^{2}}(v^{2}-1)$ and $d=\frac{\pi^{2}}{8m_{0}v}$. In
Figs. 5(a,b), we plot the angle of rotation of bases
$\phi$($=\frac{\Psi}{2}$), using the perturbed soliton given in Eq.(30) for
the same parametric choices as in the case of localized inhomogeneity. We
observe from the figures that, the fluctuation appears in the width of the
soliton without any change in its topological character asymptotically. In
order to compare the results with that of the case when the hydrogen bonding
inhomogeneity is absent as also found in [28](see Figs. (6a,b)), we plot the
perturbed soliton found in (30) when $B=0$. As in the previous case, here also
we observe that fluctuation appears in the width of the soliton in both the
cases and the inhomogeneity in hydrogen bonds introduces more fluctuation.

##  4 Conclusions

In this paper, we studied the nonlinear dynamics of a completely inhomogeneous
(inhomogeneity in both stacking and hydrogen bonding) DNA double helix by
considering the dynamical plane-base rotator model. The dynamics of this model
in the continuum limit gives rise to a perturbed sine-Gordon equation, which
was derived from the Hamiltonian consisting of site-dependent stacking and
hydrogen bonding energies. In the unperturbed limit, the dynamics is governed
by the kink-antikink soliton of the integrable sine-Gordon equation which
represents the opening of base pairs in a homogeneous DNA. In order to
understand the effect of inhomogeneity in stacking and hydrogen bonds on the
base pair opening, we carried out a perturbation analysis using multiple-scale
soliton perturbation theory. The perturbation not only modifies the shape of
the soliton but also introduces change in the velocity of the soliton. From
the results, we observe that when the inhomogeneity is in a localized or
periodic form, the width of the soliton remains constant. However, the
velocity of the soliton increases, decreases or remains uniform and even the
soliton stops, depending on the values of inhomogeneity represented by $A$ and
$B$. The soliton in all the cases are found to be stable. From the results of
the perturbed soliton we observe that the inhomogeneity in stacking and
hydrogen bonds in both the cases (localized and periodic forms) introduce
fluctuation in the form of pulses in the width of the soliton. However, there
is no change in the topological character of the soliton in the asymptotic
region(see Figs. 3 and 5). The fluctuation is more when both the
inhomogeneities (site-dependent stacking and hydrogen bonds) are present,
whereas it is less in the case of homogeneous hydrogen bonds and site-
dependent stacking. Hence, we conclude that, the addition of inhomogeneity in
hydrogen bonds does not introduce big changes in the soliton parameters and
shape except a correction in the velocity of the soliton and fluctuation. The
above dynamical behaviour may act as energetic activators of the RNA-
polymerase transport process during transcription in DNA. In the case of short
DNA chains the discreteness effect assumes importance and hence we will
analyse the discrete dynamical equations (2) and the results will be published
elsewhere.

##  5 Acknowledgments

The work of M. D forms part of a major DST project. V. V thanks SBI for
financial support.

## References

  * [1] L. V. Yakushevich, Physica D 79 (1994) 77.
  * [2] L. V. Yakushevich, Nonlinear Physics of DNA ( Wiley-VCH, Berlin, 2004).
  * [3] S. W. Englander, N. R. Kallenbanch, A. J. Heeger, J. A. Krumhansl and S. Litwin, Proc. Natl. Acad. Sci. U.S.A 77 (1980) 7222.
  * [4] S. Yomosa, Phys. Rev. A 27 (1983) 2120.
  * [5] S. Yomosa, Phys. Rev. A 30 (1984) 474.
  * [6] S. Takeno and S. Homma, Prog. Theor. Phys. 70 (1983) 308.
  * [7] S. Takeno and S. Homma, Prog. Theor. Phys. 72 (1984) 679.
  * [8] S. Takeno, Phys. Letts. A 339 (2005) 352.
  * [9] P. Jensen, M. V. Jaric and K. H. Bannenmann, Phys. Letts. A 95 (1983) 204.
  * [10] A. Khan, D. Bhaumik and B. Dutta-Roy, Bull. Math. Biol. 47 (1985) 783.
  * [11] V. K. Fedyanin and V. Lisy, Studia Biophys. 116 (1986) 65.
  * [12] R. V. Polozov and L. V. Yakushevich, J. Theor. Biol. 130 (1988) 423.
  * [13] J. A. Gonzalez and M. M. Landrove, Phys. Letts. A 292 (2002) 256.
  * [14] L. V. Yakushevich, Nanobiology 1 (1992) 343 .
  * [15] M. Peyrard and A. R. Bishop, Phys. Rev. Lett. 62 (1989) 2755.
  * [16] G. Kalosakas, K. Q. Rasmussen and A. R. Bishop, Synthetic Metals, 141 (2004) 93.
  * [17] P. L. Christiansen, P. C. Lomdahl and V. Muto, Nonlinearity 4 (1991) 477.
  * [18] T. Dauxios , Phys. Lett. A, 159 (1991) 390.
  * [19] C. B. Tabi, A. Mohamdou and T. C. Kofane, Phys. Scr. 77 (2008) 045002.
  * [20] M. Cadoni and R. De Leo, Phys . Rev. E, 021919 (2007).
  * [21] M. Barbi, S. Cocco and M. Peyrard, Phys. Letts. A, 253 (1999) 358.
  * [22] M. Barbi, S. Cocco , M. Peyrard and N. Theodorakopoulos, Phys. Rev. E, 68 (2003) 061909.
  * [23] S. Cocco and R. Manasson, Phys. Rev. Lett. 83 (1999) 5178.
  * [24] A. Campa, Phys. Rev. E , 63 (2001) 021901.
  * [25] J. Ladik and J. Cizek, Int. J. Quantum Chem. 26 (1984) 955.
  * [26] E. Cubero, E. C. Sherer, F. J. Luque, M. Orozco and C. A. Laughton, J. Am. Chem. Soc. 121 (1999) 8653 .
  * [27] M. Hisakado and M. Wadati, J. Phys. Soc. Jpn. 64 (1995) 1098.
  * [28] M. Daniel and V. Vasumathi, Physica D 231 (2007) 10.
  * [29] M. J. Ablowitz, D. J. Kaup, A. C. Newell, and H. Segur, Stud. Appl. Math. 53 (1974) 249.
  * [30] J. Yan and Y. Tang, Phys. Rev. E 54 (1996) 6816.
  * [31] Y. Tang and W. Wang, Phys. Rev. E 62 (2000) 8842.
  * [32] J. Yan, Y. Tang, G. Zhou and Z. Chen, Phys. Rev. E 58 (1998) 1064.
  * [33] Y. Matsuno, J. Math. Phys. 41 (2000) 7061.
  * [34] F. Zhang, Y. S. Kivshar and L. Vazquez, Phys. Rev. A 45 (1992) 6019.
  * [35] L. V. Yakushevich, A. V. Savin and L. I. Manevitch, Phys. Rev. E 66 (2002) 016614.
  * [36] R. Dandoloff and A. Saxena, J. Phys.: Condens. Matter 9 (1997) L667.
  * [37] A. Sanchez, A. R. Bishop and F. Dominguez-Adame, Phys. Rev. E 49 (1994) 4603.
  * [38] E. Kreyszig, Advanced Engineering Mathematics (John-Wiley, New York, 2002).

Fig.1. (a) A schematic structure of the B-form DNA double helix. (b) A
horizontal projection of the $n^{th}$ base pair in the xy-plane.
Fig.2. (a) Kink and (b) antikink soliton solutions of the sine-Gordon equation
(Eq.(4) when $\epsilon=0$ ) with $v_{0}=0.4$. (c) A sketch of the formation of
open state configuration in terms of kink-antikink soliton in DNA double
helix.
Fig.3.(a) The perturbed kink-soliton and (b) the perturbed antikink-soliton
for the inhomogeneity $g(z)=h(z)=sechz$ with $A=B=1$ and $v_{0}=0.4$.
Fig.4. (a) The perturbed kink-soliton and (b) the perturbed antikink-soliton
for the inhomogeneity $s(z)=sechz,h(z)=0$ with $A=1$ and $v_{0}=0.4$.
Fig.5.(a) The perturbed kink-soliton and (b) the perturbed antikink-soliton
for the inhomogeneity $g(z)=h(z)=\cos z$ with $A=B=1$ and $v_{0}=0.4$.
Fig.6.(a) The perturbed kink-soliton and (b) the perturbed antikink-soliton
for the inhomogeneity $g(z)=\cos z,h(z)=0$ with $A=1$ and $v_{0}=0.4$.

Figure 1: (a) A schematic structure of the B-form DNA double helix. (b) A
horizontal projection of the $n^{th}$ base pair in the xy-plane. Figure 2: (a)
Kink and (b) antikink soliton solutions of the sine-Gordon equation (Eq.(4)
when $\epsilon=0$ ) with $v_{0}=0.4$. (c) A sketch of the formation of open
state configuration in terms of kink-antikink soliton in DNA double helix.
Figure 3: (a) The perturbed kink-soliton and (b) the perturbed antikink-
soliton for the inhomogeneity $g(z)=h(z)=sechz$ with $A=B=1$ and $v_{0}=0.4$.
Figure 4: (a) The perturbed kink-soliton and (b) the perturbed antikink-
soliton for the inhomogeneity $s(z)=sechz,h(z)=0$ with $A=1$ and $v_{0}=0.4$.
Figure 5: (a) The perturbed kink-soliton and (b) the perturbed antikink-
soliton for the inhomogeneity $g(z)=h(z)=\cos z$ with $A=B=1$ and $v_{0}=0.4$.
Figure 6: (a) The perturbed kink-soliton and (b) the perturbed antikink-
soliton for the inhomogeneity $g(z)=\cos z,h(z)=0$ with $A=1$ and $v_{0}=0.4$.


# Tactical cooperation of defectors in a multi-stage public goods game

Attila Szolnoki  Xiaojie Chen  Institute of Technical Physics and Materials
Science, Centre for Energy Research, P.O. Box 49, H-1525 Budapest, Hungary
School of Mathematical Sciences, University of Electronic Science and
Technology of China, Chengdu 611731, China

###### Abstract

The basic social dilemma is frequently captured by a public goods game where
participants decide simultaneously whether to support a common pool or not and
after the enhanced contributions are distributed uniformly among all
competitors. What if the result of common efforts is not distributed
immediately, but it is reinvested and added to the pool for a next round? This
extension may not only result in an enhanced benefit for group members but
also opens new strategies for involved players because they may act in
distinct rounds differently. In this work we focus on the simplest case when
two rounds are considered, but the applied multiplication factors dedicated to
a certain round can be different. We show that in structured populations the
winning strategy may depend sensitively on the ratio of these factors and the
last round has a special importance to reach a fully cooperative state. We
also observe that it may pay for defectors to support the first round and
after enjoy the extra benefit of accumulated contributions. Full cooperator
strategy is only viable if the second round ensures a premium benefit of
investments.

###### keywords:

public goods game , cooperation , phase transition

††journal: Chaos Solitons and Fractals

##  1 Introduction

According to the evolutionary selection principle, when individual interests
are in conflict the most successful strategy prevails among competing
approaches [1]. For an adequate mathematical description we introduce a
measure of success, which is a payoff value gained by a player when
interacting with others. In other words, the collected payoff is a clear
feedback for a player to decide on which strategy to choose. Importantly, this
assumption still allows several ways how to implement the above mentioned
selection principle because many alternative microscopic dynamical rules can
serve this goal [2, 3, 4]. Beside the frequently applied imitating the more
successful partner rule, we may apply alternative standards of social
learning, like using myopic update, birth-death, or death-birth process, etc.
[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. For further details our reader is
advised to check related topical reviews [16, 17, 18].

In the standard public goods game ($PGG$), which is a simple metaphor of the
conflict of two major strategies, a player decides whether to contribute to a
common pool or not [19, 20]. When decisions are made, we enhance the
accumulated contributions by a multiplication factor and redistribute it among
all group members uniformly independently of their strategies [21, 22, 23].
Therefore, it is not surprising that those who contribute first gain less than
the defector players who just enjoy the fruit of others’ efforts. This
annoying conclusion cannot be avoided, just only if we assume some additional
circumstance, like punishing defectors, rewarding cooperators [24, 25, 26, 27,
28, 29], or other more sophisticated mechanism which allows cooperators to
collect competitive payoff [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40].

From our present viewpoint, however, it is more important when to realize the
consequence of a chosen strategy, because players may not necessarily face
with the result of their choices immediately. For example, it can happen that
there is some delay when the results of strategy interactions become available
[41, 42]. Seasonal effects or long-term investments may explain why we not
always evaluate strategy success immediately [43]. Alternatively, there are
situations when the fruit of collective efforts is not distributed at once,
but it is considered as an additional contribution to a new round when players
are invited to cooperate to the common pool again. In this paper we study the
possible consequences of such multi-stage $PGG$ by analyzing a two-level
version of the classic game. This extension makes players possible to choose
alternative strategies in different rounds, which enlarges the possible number
of competing strategies. As an important extension, we do not insist that the
enhancement factors dedicated to alternating rounds should be equal, but they
could be different. This liberty helps us to reveal that the ratio of the
mentioned multiplications factors has a special importance if players are
arranged in a spatially structured population. Furthermore, the success of the
application of mixed strategies is highly biased, because in certain parameter
region it could be useful for defectors to contribute to the first round and
enjoy the accumulated benefits in the last round.

In the next section we specify the details of the extended $PGG$ and discuss
the evolutionary outcome in a mixed population where all players can interact
with anybody else randomly. Our key observations are summarized in Sec. 3
where we present the phase diagram of a spatially structured population and
give further details of the emerging phase transitions. In the last section we
briefly analyze the possible general consequences of multi-stage public goods
game and briefly discuss further potential extensions.

##  2 Multi-stage public goods game

By following the standard definition of $PGG$ we consider a population of $N$
players who can be cooperators ($C$) or defectors ($D$). $G$ number of players
form a group to execute a collective project. While cooperators contribute a
$c=1$ amount to the common pool, defectors avoid such effort. Instead, they
only enjoy the benefit of the joint venture. The sum of collective investment
is multiplied by an $r$ enhancement factor which expresses the fact that
collective efforts result in higher income level than the simple sum of
individual contributions. This enhanced amount is distributed among all group
members equally, no matter if a player properly contributed to the joint
effort or not. Regarding the microscopic dynamics, we apply the widely
accepted pairwise comparison imitation rule. In particular, a player $y$ who
has strategy $s_{y}$ adopts the $s_{x}$ strategy of a player $x$ with the
following $\Gamma$ probability:

$\Gamma(s_{x}\to s_{y})=[1+\exp(\Pi_{s_{y}}-\Pi_{s_{x}})/K]^{-1}\,.$

Here $\Pi_{s_{y}}$ and $\Pi_{s_{x}}$ denote the payoff values of involved
strategies, while $K$ parameter quantifies the amplitude of noise level during
the adoption process. For proper comparison with previous studies of standard
model we here apply $K=0.5$ value which ensures a likely adoption of more
successful strategy, but also allows a reverse process with a small
probability [44].

Next we extend the standard model and introduce two-stages of the game to
calculate the proper payoff values of interacting players. More precisely,
when the involved group members contribute to the common pool and their
investments are enhanced by an $r_{1}$ multiplication factor then we do not
distribute this sum. Instead, we invest the whole amount into a second pool to
increase the additional contributions of group members. To distinguish the
potential weights of various stages, in the second round we apply $r_{2}$
multiplication factor to enhance the sum of players’ contributions and the
amount saved from the first round. Importantly, a player may decide to
contribute to the first or/and second round, hence we can distinguish four
different strategies. They are designated as $DD$, $CD$, $DC$, and $CC$
depending on whether a player contributes to the first or/and to the second
round, or not. Accordingly, the payoff value collected by different strategies
from a group interaction are the following:

$\displaystyle\Pi_{DD}$ $\displaystyle=$
$\displaystyle\frac{r_{2}(r_{1}[n_{CD}+n_{CC}]+[n_{DC}+n_{CC}])}{G}\,,$
$\displaystyle\Pi_{CD}$ $\displaystyle=$ $\displaystyle\Pi_{DD}-1\,,$
$\displaystyle\Pi_{DC}$ $\displaystyle=$ $\displaystyle\Pi_{DD}-1\,,$
$\displaystyle\Pi_{CC}$ $\displaystyle=$ $\displaystyle\Pi_{DD}-2\,,$

where $n_{s}$ denotes the number of players having strategy $s$ in the group.
Here $DD$ can be considered as a “classic” defector who never contributes to
joint ventures, while a $CC$ player invests to both pools, hence can be
considered as an unconditional cooperator. A key question is whether it pays
to mix the attitudes by cooperating (or defect) in the first round and apply
the reversed act in the second round. Or, is there a coexistence of certain
strategies at specific parameter values?

We consider the generalized problem both in well-mixed and in structured
populations. While the former is mathematically feasible in several cases, but
the latter option in general is significantly closer to realistic situations
[45, 46, 47, 48, 49, 50, 51] Therefore it is always instructive to compare the
system behaviors for both conditions [52, 53, 54]. For a proper comparison
with previous results of traditional $PGG$ we here apply a square lattice
topology to describe a spatially structured population where every player has
four nearest neighbors, hence they form a $G=5$ member group [44]. Of course,
a player is involved not only in one game, but also in four others ones which
are organized by the mentioned neighbors. To gain reliable statistic we used
$L\times L$ = $800\times 800$ system size where the requested relaxation time
to reach the stationary state varied between $10^{4}$ to $10^{6}$ Monte Carlo
steps ($MCS$), depending on the proximity of phase transition points.
According to the standard simulation protocol, during a full $MCS$ every
player has a chance on average to adopt a strategy from a partner.

Before presenting the subtle system behavior of spatially structured
population we first summarize the result for a well-mixed population. To make
the comparison accurate we here assume that a player is also involved in $G$
games as for the spatial system. Perhaps it is worth noting that in the
traditional $PGG$ defectors prevail if $r$ is lower than $G$ and cooperators
can survive only above this threshold value. In our two-stage game the system
behavior is conceptually similar, as shown in Fig. 1. Here we plotted the
evolutionary stable solutions in dependence of the multiplication factors
dedicated to the specific stage of the game. Similarly to the traditional one-
stage game persistent defectors, i.e. $DD$ strategies are always selected if
the product of $r_{1}\cdot r_{2}$ does not exceed $G\cdot G$ threshold value.
Above this point $CD$ strategy prevails.

Figure 1: Phase diagram of a mixed population on the plane of multiplication
factors which are dedicated to specific rounds of $PGG$. If these factors are
low, the system evolves into a trapped phase where only $DD$ players are
present. If the product of $r_{1}\cdot r_{2}$ exceeds $G\cdot G$ then $CD$
strategy prevails. These phases are separated by a discontinuous phase
transition, marked by dashed blue line.

In other words, when the product of enhancement factors is high enough then it
pays for defector players to invest into the first round and after enjoy the
accumulated benefit of the two-stage game. It is also worth noting that the
liberty to choose diverse values of enhancement factor for different stages
has no any relevance: the roles of enhancement factors are symmetric and the
new solution can be reached independently of which multiplication factor is
higher.

##  3 Phase transitions in a spatially structured population

In this section we give a detailed report about how a structured population
behaves when a two-stage $PGG$ is considered. In evolutionary game theory we
frequently observe that the system behavior is strikingly different when we
leave the analytically feasible well-mixed condition and assume an interaction
topology where players have fixed and limited number of partners [55, 56, 57].
This is the case in our present model, too, as it is demonstrated in Fig. 2.

Figure 2: Phase diagram of the system behavior in case the players are
arranged on a square lattice. Solid red lines mark continuous phase
transitions, while dashed blue lines denote discontinuous phase transition
points. If we compare this diagram with the one obtained for well-mixed
population then we can observe a significantly richer behavior where seven
evolutionary stable solutions can be found for different parameter pairs.

The first observation is the large diversity of stable solutions on the
parameter plane. Namely, we can detect seven different phases as we change the
($r_{1},r_{2}$) parameter pair. Similarly to the well-mixed case $DD$ strategy
becomes dominant for small values of enhancement factors, but this phase is
significantly smaller here. The tactical cooperation of defectors, in other
words the $CD$ strategy, can also be the winner of the evolutionary process,
but only if $r_{2}$ value remains moderate. Furthermore, persistent
cooperation of $CC$ strategy could be dominant for high $r_{2}$ values
independently of the value of $r_{1}$. This feature underlines the strong
asymmetry in the role of enhancement factors, which cannot be detected in a
well-mixed population.

Beside the mentioned pure phases we can observe coexistence of different
strategies. For example $DD$ and $CD$ strategies can form a stable solution
between the related pure phases. Similar behavior can be observed for the
relation of $CD$ and $CC$ strategies or the coexistence of the extreme
strategies of $DD$ and $CC$. Interestingly, $DD$, $CD$, and $CC$ strategies
can form a stable three-member solution at specific pairs of enhancement
factors. But $DC$ strategy is not proved to be viable for any combination of
parameters.

The mentioned phases are separated by continuous (marked by solid red) or
discontinuous phase transition points. The latter borders are denoted by
dashed blue lines. To illustrate the richness of system behavior in Fig. 3 we
present a cross-section of the diagram obtained at $r_{1}=1.16$ as we
gradually increased the value of $r_{2}$. By solely varying $r_{2}$ we can
detect six consecutive phase transitions. Their positions are marked by arrows
on the top. The first transition happens at $r_{2}=3.226$ when $CD$ becomes
viable and coexist with $DD$ strategy due to the relatively high value of
$r_{2}$. After, at $r_{2}=3.415$ an interesting first-order phase transition
happens, where a qualitatively new solution, the coexistence of $DD$ and $CC$
emerges. Later we will give further details which explains this discontinuous
phase transition.

As we increase $r_{2}$ further, $CD$ strategy becomes viable again and forms a
stable three-member solution with the mentioned strategies from $r_{2}=4.488$.
However, as the phase diagram highlights, this solution is restricted to a
very narrow area of parameters because it requires a not too large $r_{1}$
value, otherwise $CD$ would enjoy its first investment too easily. In
parallel, $r_{2}$ should be at an intermediate value: for smaller $r_{2}$
values $DD$ would dominate, while a high second enhancement factor would help
cooperator strategies. This is why at $r_{2}=4.676$ cooperator strategies,
$CD$ and $CC$ become more efficient and form a more successful two-strategy
solution again. Lastly, we can observe a rather exotic reentrant to $DD+CC$
phase, but this recovery of the mentioned solution is strongly related to the
low value of $r_{1}$. For larger $r_{1}$ values, when the first round is more
effective, this solution cannot be detected anymore, but gives way to the $CD$
tactical defector strategy. Staying at the mentioned exotic phase transitions,
the first one is discontinuous at $r_{2}=5.404$ while the second one is
continuous at $r_{2}=5.535$, as it is illustrated in the inset of Fig. 3 where
we zoomed the portion of $DD$ strategy for clarity.

Figure 3: The portion of surviving strategies in dependence of $r_{2}$
enhancement factor as we cross the phase diagram at a fixed $r_{1}=1.16$
value. The gradual change of $r_{2}$ results in six(!) consecutive phase
transitions whose positions are marked by arrows on the top. The last two
transition points are hardly visible on this scale therefore we zoomed the
portion of $DD$ strategy in the inset. The error bars are comparable to the
thickness of the lines.

Of course, we can detect a series of phase transitions as we increase the
benefit of first round at a fixed $r_{2}$ value, but the system behavior is
simpler than previously. A representative cross-section is given in Fig. 4
where we plot the portion of surviving strategies as a function of $r_{1}$ at
a fixed $r_{2}=3.7$ value. In the low $r_{1}$ region only the extreme
strategies, $DD$ and $CC$, can survive. This fact can be interpreted that at
very low $r_{1}$ the system becomes practically identical to a classical,
single-stage $PGG$, where only pure strategies can be viable. Their relations
depends only on the value of $r_{2}$, which takes the role of the traditional
enhancement factor, $r$. When it is low then only $DD$ survives, while high
$r_{2}$ can provide the full domination of $CC$. Between these extreme cases
the mentioned strategies coexist due to the spatiality of the population. This
is the case at $r_{2}=3.7$, as we mentioned above. One may note that this
value is below the $r=3.745$ threshold value of classic model [44], but of
course here the first round may contribute a bit which lowers this critical
point modestly. Turning back to Fig. 4, as we increase $r_{1}$, $CD$ becomes
viable and form the previously mentioned three-member solution with their
partners. By supporting the first round further, $CD$ becomes really effective
and crowds out classic $CC$ strategy. In other words, it pays defectors to
invest into the first round because later they can harvest a significantly
higher payoff at the end. Figure 4 also shows that this tactical defector
strategy becomes exclusive at a relatively small $r_{1}$ value by giving no
space even for $DD$ strategy. The phase diagram shown in Fig. 2 also
highlights that this strategy remains viable for the majority of
($r_{1},r_{2}$) parameter pairs and only very low $r_{1}$ or very high $r_{2}$
is capable to replace it with a traditional strategy.

Figure 4: The portion of surviving strategies in dependence of $r_{1}$, as we
cross the phase diagram at fixed $r_{2}=3.7$ value. Three phase transitions
can be detected as we reach the full $CD$ phase at high $r_{1}$ values. Their
positions are marked by arrows. The error bars are comparable to the thickness
of the lines.

In the following we give a deeper insight about the origin of discontinuous
phase transitions we previously reported. For this reason we monitor the
strategy evolution at two parameter pairs which are very close to each other,
but they are in different sides of the transition point. More precisely, we
follow the evolution at $r_{1}=1.16$, $r_{2}=3.41$ and at $r_{1}=1.16$,
$r_{2}=3.42$. To make our point more visible, we use a prepared initial state
where we divide the available space into two halves where the first half
contains $DD$ and $CC$ players randomly, while the second half is occupied by
$DD$ and $CD$ players. This initial setup is illustrated in panel (a) of Fig.
5.

Here, in agreement with the color code used in Fig. 3 and in Fig. 4, we marked
by red color $DD$ players, while blue (green) color denotes $CC$ ($CD$)
players on the grid. After we launch the evolution according to the two-stage
$PGG$, but imitations are allowed only within the subsystems we defined. As a
result, in both halves of the space a stationary solution emerges. The
mentioned stage is illustrated in Fig. 5(b). This panel highlights that both
$DD+CC$ and $DD+CD$ could be a destination of the evolutionary process. We
note that the presented patterns were obtained at $r_{2}=3.41$, but
practically we would get identical strategy distribution for $r_{2}=3.42$ by
following similar protocol. These solutions can be reached after 5000 $MCS$ of
relaxation. In the following we open the borders for imitations and allow the
solutions to compete for space. Accordingly, either $DD+CC$ or $DD+CD$ will
dominate and invade the whole space depending on the actual value of $r_{2}$.
These final destinations are not shown here but they show similar patterns as
the sub-system solutions depicted in panel (b).

Figure 5: Stability analysis of sub-solutions to explain the origin of
discontinuous phase transitions obtained at $r_{2}=3.41$ (c) and $r_{2}=3.42$
(d) in a $800\times 800$ system. In both cases $r_{1}=1.16$. Further details
and more specific explanation can be found in the main text.

The trajectories of the reported evolutions are shown in panel (c) and panel
(d), where $r_{2}=3.41$ and $r_{2}=3.42$ values were used respectively. Here
the random initial states of the subsystems were started at -5000 $MCS$, and
the borders were opened at zero time, which are marked by arrows in these
panels. As panel (c) illustrates, the portion of $CC$ strategy starts falling
immediately, while it grows when we gently modified the $r_{2}$ value in panel
(d). Finally the system evolves into the mentioned two-member solutions when
the third strategy dies out. We note that we used $800\times 800$ system size
to avoid too noisy trajectories, but panel (a) and (b) show only just a
$200\times 200$ portions around the borders, otherwise the compact domains of
$CC$ and $CD$ strategies would be less visible.

This comparison underlines that both solutions would be a stable outcome in
the absence of the other one. Technically, both $CC$ and $CD$ players form
compact islands in the sea of $DD$ players, as it can be seen in panel (b).
This is how network reciprocity mechanism works in general. Namely, these
convex patterns make the cooperating strategies possible to reach competitive
payoff despite of the relatively low values of enhancement factors. When we
are at the phase transition point exactly then these solutions are in
equilibrium. But if we leave this point slightly then this equilibrium is
broken and one of the solutions dominates the other. Notably, in the vicinity
of this transition point the relaxation is slow, as it is illustrated in panel
(c) and (d) in Fig. 5, because the driving force to reach the final state is
weak. But this destination is inevitable if we wait enough in a necessary
large system. Summing up, what we can observe here are the typical features
which characterize first-order phase transitions in statistical physics.

Last we turn back to the comparison of outcomes obtained in well-mixed and in
structured populations. As we already noted, a significant difference can be
detected in the roles of enhancement factors. While in a random system it has
no particular significance which factor is larger because their product
counts, in structured population this symmetry is seriously broken. To support
our argument in Fig. 6 we show the average cooperation level of players during
the two-stage game. In agreement with previously used color coding we here
mark by red color those parameter pairs where players never cooperate in the
stationary state. Similarly, blue color denotes those parameter values when
players behave as classical cooperators and contribute to both pools. Green
color shows those $(r_{1},r_{2})$ pairs when there is one contribution on
average. Perhaps it is worth noting that this average can be reached in
different ways: a mixture of $DD$ and $CC$ players or pure $CD$ players can
provide similar levels of cooperation. The main point here is the strong
asymmetry of heat map in terms of the values of enhancement factors. In
particular, the value of $r_{1}$ could be really high, but if $r_{2}$ is
modest, we can only reach an intermediate cooperation level. On the other
hand, if $r_{2}$ is high enough then we can reach the desired high cooperation
level almost independently of the value of $r_{1}$.

Figure 6: Heat map of cooperation level on the parameter plane. Red color
denotes the fully defector state while blue color marks full cooperation.
Green color shows an intermediate state when players partly cooperate on
average during the two-stage process. Legend on the right quantifies these
levels. The asymmetry how the final outcome depends on the values of
enhancement factors is striking.

##  4 Conclusions

It is always an individual decision whether we want to cooperate and
contribute to a common venture or not. In short term it could pay to save such
cost and only enjoy the fruit of efforts, but in long term, if everybody else
think similarly then we cannot avoid the so-called tragedy of the commons
trajectory [58, 59, 60, 61, 62, 63, 64]. In other words, we will face to the
consequence of our choice. But this outcome should not be clear immediately if
we are part of a multi-stage project and the success of our choice, namely the
collected payoff will only be clear at the very end. A delayed evaluation or
the consequence of a multi-level project, however, opens a new direction for
players to apply more sophisticated strategies than the classic unconditional
cooperator and defector choices. For example, they may cooperate or defect in
the early stage and later they may alter their attitudes. In this work we
explored the consequence of a two-stage public goods game, where the fruit of
the first round is not distributed among the players, but it is reinvested
into the second stage. It is an important extension that the weight of
different rounds could be unequal, which also opens to study the combination
of heterogeneous joint ventures.

We have shown that the proposed extension has no particularly interesting
consequence in a well-mixed population because it simply rescale the
enhancement factors. Therefore only unconditional defectors survive if the
product of enhancement factors remains below a threshold value. The only
interesting feature is emergence of tactical defectors who dominate above the
mentioned threshold value. They only invest onto the first round, and enjoy
the accumulated benefit. The roles of enhancement factors dedicated to a
certain round remain symmetric and only their product counts.

This simple picture becomes more complex in structured populations. We here
presented the results obtained on square grid which makes also possible to
compare the results with traditional model played on similar topology. The
presented phase diagram revealed several phases which are formed by a pure
strategy or the combination of these. We could detect not only two-member
coexistence, but also a solution in which three strategies can survive in the
stationary state. These phases are separated by continuous or discontinuous
phase transition lines. We have given explanation why the latter can emerge in
this system.

Perhaps the most spectacular feature of the spatially structured population is
the strong symmetry breaking regarding the roles of enhancement factors
dedicated to certain rounds. A moderate level of cooperation can be reached at
relatively low $r_{2}$ value if the $r_{1}$ value is high enough, but complete
cooperation can only be reached for a high $r_{2}$ value. In the latter case,
however, the value of $r_{1}$ has no particular importance. This observation
may help us how to design multi-stage games because the last round has
paramount importance to reach a high cooperation level.

To follow this path we may study what if more stages are present in a $PGG$
where the results of previous rounds are reinvested again and again until the
last stage when all cards are revealed. Of course, such an extension makes the
study really challenging because the number of the available strategies grows
exponentially with the number of rounds. Furthermore, if we allow that
distinct stages may influence the final outcome in different ways, i.e. we
introduce stage-specific enhancement factors, then the optimal player’s
strategy could be more interesting. For example, successful players may want
to defect, cooperate and defect again. Hopefully, following works will clarify
this research avenue.

This research was supported the National Natural Science Foundation of China
(Grants No. 61976048 and No. 62036002).

## References

  * Maynard Smith [1982] J. Maynard Smith, Evolution and the Theory of Games, Cambridge University Press, Cambridge, U.K., 1982.
  * Amaral and Javarone [2020] M. A. Amaral, M. A. Javarone, Strategy equilibrium in dilemma games with off-diagonal payoff perturbations, Phys. Rev. E 101 (2020) 062309.
  * Yang and Yang [2019] H.-X. Yang, J. Yang, Reputation-based investment strategy promotes cooperation in public goods games, Physica A 523 (2019) 886–893.
  * Jiao et al. [2020] Y. Jiao, T. Chen, Q. Chen, The impact of expressing willingness to cooperate on cooperation in public goods game, Chaos, Solit. Fract. 140 (2020) 110258.
  * Szabó and Szolnoki [2012] G. Szabó, A. Szolnoki, Selfishness, fraternity, and other-regarding preference in spatial evolutionary games, J. Theor. Biol. 299 (2012) 81–87.
  * Liu et al. [2019a] R.-R. Liu, C.-X. Jia, Z. Rong, Effects of enhancement level on evolutionary public goods game with payoff aspirations, Appl. Math. Comput. 350 (2019a) 242–248.
  * Szolnoki and Chen [2018] A. Szolnoki, X. Chen, Competition and partnership between conformity and payoff-based imitations in social dilemmas, New J. Phys. 20 (2018) 093008.
  * Cardoso et al. [2020] F. M. Cardoso, C. Gracia-Lázaro, Y. Moreno, Dynamics of heuristics selection for cooperative behaviour, New J. Phys. 22 (2020) 123037.
  * Zhang et al. [2019] L. Zhang, C. Huang, H. Li, Q. Dai, Aspiration-dependent strategy persistence promotes cooperation in spatial prisoner’s dilemma game, EPL 126 (2019) 18001.
  * Fu et al. [2019] M. Fu, W. Guo, L. Cheng, S. Huang, D. Chen, History loyalty-based reward promotes cooperation in the spatial public goods game, Physica A 525 (2019) 1323–1329.
  * Szolnoki and Danku [2018] A. Szolnoki, Z. Danku, Dynamic-sensitive cooperation in the presence of multiple strategy updating rules, Physica A 511 (2018) 371–377.
  * Xu et al. [2019] Z. Xu, R. Li, L. Zhang, The role of memory in human strategy updating in optional public goods game, Chaos 29 (2019) 043128.
  * Li et al. [2021a] X. Li, L. Cheng, X. Niu, S. Li, C. Liu, P. Zhu, Highly cooperative individuals’ clustering property in myopic strategy groups, Eur. Phys. J. B 94 (2021a) 126\.
  * Lin et al. [2020] J. Lin, C. Huang, Q. Dai, J. Yang, Evolutionary game dynamics of combining the payoff-driven and conformity-driven update rules, Chaos, Solitons and Fractals 140 (2020) 110146.
  * Zhang et al. [2021] L. Zhang, C. Huang, H. Li, Q. Dai, J. Yang, Cooperation guided by imitation, aspiration and conformity-driven dynamics in evolutionary games, Physica A 561 (2021) 125260.
  * Szabó and Fáth [2007] G. Szabó, G. Fáth, Evolutionary games on graphs, Phys. Rep. 446 (2007) 97–216.
  * Perc et al. [2013] M. Perc, J. Gómez-Gardeñes, A. Szolnoki, L. M. Floría and Y. Moreno, Evolutionary dynamics of group interactions on structured populations: a review, J. R. Soc. Interface 10 (2013) 20120997.
  * Perc et al. [2017] M. Perc, J. J. Jordan, D. G. Rand, Z. Wang, S. Boccaletti, A. Szolnoki, Statistical physics of human cooperation, Phys. Rep. 687 (2017) 1–51.
  * Sigmund [2010] K. Sigmund, The Calculus of Selfishness, Princeton University Press, Princeton, NJ, 2010.
  * Chen and Szolnoki [2016] X. Chen, A. Szolnoki, Individual wealth-based selection supports cooperation in spatial public goods game, Sci. Rep. 6 (2016) 32802.
  * Liu et al. [2021] J. Liu, M. Peng, Y. Peng, Y. Li, C. Chu, X. Li, Q. Liu, Effects of inequality on a spatial evolutionary public goods game, Eur. Phys. J. B 94 (2021) 167.
  * Quan et al. [2021a] J. Quan, C. Tang, X. Wang, Reputation-based discount effect in imitation on the evolution of cooperation in spatial public goods games, Physica A 563 (2021a) 125488.
  * Rong et al. [2019] Z. Rong, Z.-X. Wu, X. Li, P. Holme, G. Chen, Heterogeneous cooperative leadership structure emerging from random regular graphs, Chaos 29 (2019) 103103.
  * Chen et al. [2014] X. Chen, A. Szolnoki, M. Perc, Probabilistic sharing solves the problem of costly punishment, New J. Phys. 16 (2014) 083016.
  * Lucas S. Flores [2021] M. A. A. M. H. V. Lucas S. Flores, Heitor C.M. Fernandes, Symbiotic behaviour in the public goods game with altruistic punishment, J. Theor. Biol. 524 (2021) 110737.
  * Chen and Szolnoki [2018] X. Chen, A. Szolnoki, Punishment and inspection for governing the commons in a feedback-evolving game, PLoS Comput. Biol. 14 (2018) e1006347.
  * Lv and Song [2022] S. Lv, F. Song, Particle swarm intelligence and the evolution of cooperation in the spatial public goods game with punishment, Appl. Math. Comput. 412 (2022) 126586.
  * Fu et al. [2021a] M. Fu, J. Wang, L. Cheng, L. Chen, Promotion of cooperation with loyalty-based reward in the spatial prisoner’s dilemma game, Physica A 580 (2021a) 125672.
  * Yang and Fu [2020] H.-X. Yang, M.-J. Fu, A punishment mechanism in the spatial public goods game with continuous strategies, EPL 132 (2020) 10007.
  * Cong et al. [2017] R. Cong, Q. Zhao, K. Li, L. Wang, Individual mobility promotes punishment in evolutionary public goods games, Sci. Rep. 7 (2017) 14015.
  * Liu et al. [2018] J. Liu, T. Li, W. Wang, N. Zhao, F. Hang, Impact of strategy-neutral rewarding on the evolution of cooperative behavior, Chaos, Solit. and Fract. 106 (2018) 76–79.
  * Chen et al. [2015] X. Chen, A. Szolnoki, M. Perc, Competition and cooperation among different punishing strategies in the spatial public goods game, Phys. Rev. E 92 (2015) 012819\.
  * Quan et al. [2019] J. Quan, X. Li, X. Wang, The evolution of cooperation in spatial public goods game with conditional peer exclusion, Chaos 29 (2019) 103137.
  * Cheng et al. [2020] F. Cheng, T. Chen, Q. Chen, Rewards based on public loyalty program promote cooperation in public goods game, Appl. Math. and Comput. 378 (2020) 125180\.
  * Szolnoki and Perc [2021] A. Szolnoki, M. Perc, The self-organizing impact of averaged payoffs on the evolution of cooperation, New J. Phys. 23 (2021) 063068\.
  * Gao et al. [2020] S. Gao, J. Du, J. Liang, Evolution of cooperation under punishment, Phys. Rev. E 101 (2020) 062419\.
  * Li et al. [2021b] K. Li, Y. Mao, Z. Wei, R. Cong, Pool-rewarding in N-person snowdrift game, Chaos, Solit. and Fract. 143 (2021b) 110591\.
  * Wei et al. [2021] X. Wei, P. Xu, S. Du, G. Yan, H. Pei, Reputational preference-based payoff punishment promotes cooperation in spatial social dilemmas, Eur. Phys. J. B 94 (2021) 210.
  * Quan et al. [2021b] J. Quan, Z. Pu, X. Wang, Comparison of social exclusion and punishment in promoting cooperation: Who should play the leading role?, Chaos, Solitons and Fractals 151 (2021b) 111229.
  * Szolnoki and Chen [2020] A. Szolnoki, X. Chen, Blocking defector invasion by focusing on the most successful partner, Appl. Math. Comput. 385 (2020) 125430.
  * Szolnoki and Perc [2013] A. Szolnoki, M. Perc, Decelerated invasion and waning-moon patterns in public goods games with delayed distribution, Phys. Rev. E 87 (2013) 054801.
  * Yan et al. [2021] F. Yan, X. Chen, Z.Qiu, A. Szolnoki, Cooperator driven oscillation in a time-delayed feedback-evolving game, New J. Phys. 23 (2021) 053017.
  * Szolnoki and Perc [2019] A. Szolnoki, M. Perc, Seasonal payoff variations and the evolution of cooperation in social dilemmas, Sci. Rep. 9 (2019) 12575.
  * Szolnoki et al. [2009] A. Szolnoki, M. Perc, G. Szabó, Topology-independent impact of noise on cooperation in spatial public goods games, Phys. Rev. E 80 (2009) 056109\.
  * Yang et al. [2019] G. Yang, C. Zhu, W. Zhang, Adaptive and probabilistic strategy evolution in dynamical networks, Physica A 518 (2019) 99–110.
  * Broom et al. [2021] M. Broom, I. V. Erovenko, J. Rychtář, Modelling Evolution in Structured Populations Involving Multiplayer Interactions, Dyn. Games Appl. 11 (2021) 270–293.
  * Fu et al. [2021b] Y. Fu, Y. Zhang, Y. Guo, Y. Xie, Evolutionary dynamics of cooperation with the celebrity effect in complex networks, Chaos 31 (2021b) 013130.
  * Chen and Wu [2021] W. Chen, T. Wu, Adaptive multilayer networks resolve the cooperation dilemma induced by breaking the symmetry between interaction and learning, New J. Phys. 23 (2021) 093019.
  * Deng and Zhang [2021] Y. Deng, J. Zhang, The role of the preferred neighbor with the expected payoff on cooperation in spatial public goods game under optimal strategy selection mechanism, Physica A 584 (2021) 126363\.
  * Zhong et al. [2021] X. Zhong, Y. Fan, Z. Di, The evolution of cooperation in public goods games on signed networks, Physica A 582 (2021) 126217.
  * Zhu et al. [2021] P. Zhu, X. Hou, Y. Guo, J. Xu, J. Liu, Investigating the effects of updating rules on cooperation by incorporating interactive diversity, Eur. Phys. J. B 94 (2021) 58.
  * Szolnoki and Perc [2016] A. Szolnoki, M. Perc, Collective influence in evolutionary social dilemmas, EPL 113 (2016) 58004\.
  * Roca et al. [2009] C. P. Roca, J. A. Cuesta, A. Sánchez, Evolutionary game theory: Temporal and spatial effects beyond replicator dynamics, Phys. Life Rev. 6 (2009) 208–249.
  * Szolnoki and Chen [2015] A. Szolnoki, X. Chen, Benefits of tolerance in public goods games, Phys. Rev. E 92 (2015) 042813\.
  * Lütz et al. [2021] A. F. Lütz, M. A. Amaral, L. Wardil, Moderate immigration may promote a peak of cooperation among natives, Phys. Rev. E 104 (2021) 014304.
  * Li et al. [2020] Q. Li, Z. Zhang, K. Li, L. Chen, Z. Wei, J. Zhang, Evolutionary dynamics of traveling behavior in social networks, Physica A 545 (2020) 123664.
  * Zheng et al. [2021] L. Zheng, H. Xu, C. Tian, S. Fan, Evolutionary dynamics of information in the market: Transmission and trust, Physica A 581 (2021) 126228.
  * Hardin [1968] G. Hardin, The Tragedy of the Commons, Science 162 (1968) 1243–1248.
  * Cheng et al. [2019] F. Cheng, T. Chen, Q. Chen, Matching donations based on social capital in Internet crowdfunding can promote cooperation, Physica A 531 (2019) 121766\.
  * Szolnoki et al. [2011] A. Szolnoki, N.-G. Xie, C. Wang, M. Perc, Imitating emotions instead of strategies in spatial games elevates social welfare, EPL 96 (2011) 38002.
  * Liu et al. [2019b] D. Liu, C. Huang, Q. Dai, H. Li, Positive correlation between strategy persistence and teaching ability promotes cooperation in evolutionary prisoner’s dilemma games, Physica A 520 (2019b) 267–274.
  * Tao et al. [2021] Y. Tao, K. Hu, L. Shi, Risk-preference-driven participate willingness provides alternative routes to solve social dilemma, EPL 135 (2021) 28001.
  * Szolnoki et al. [2010] A. Szolnoki, Z. Wang, J. Wang, X. Zhu, Dynamically generated cyclic dominance in spatial prisoner’s dilemma games, Phys. Rev. E 82 (2010) 036110.
  * Yang and Zhang [2021] L. Yang, L. Zhang, Environmental feedback in spatial public goods game, Chaos Solit. and Fract. 142 (2021) 110485.