text:
# Dynamical mean-field approximation to coupled active rotator networks subject to white noises

## I INTRODUCTION

It has been shown that noises play intrigue and essential roles in non-linear systems. Some examples are the noise-induced transition  and the stochastic resonance . Effects of noises on the dynamics of a variety of real systems such as the Josephson junction , chemical reaction , charge density of states  and neural networks , have been studied with the use of the coupled phase model and its variants. These model are described by stochastic, nonlinear differential equations (DEs), which have been solved by simulations or analytical methods such as the Fokker-Planck equation (FPE) and the moment method. In order to make our discussion concrete, we hereafter pay our attention to the active rotator (AR) model, which was first studied by Shinomoto and Kuramoto . They obtained the FPE for the infinite-dimensional ($`N=\mathrm{}`$) AR model, discussing the phase transition between the stationary state and the time-periodic state. Responses of excitable AR models to subthreshold inputs have been investigated -. In recent years effects of additive and/or multiplicative noises in AR models have been studied by using FPE -. Although FPE method is a powerful tool for solving stochastic DEs, its application is limited to a single or infinite system subject to white noises. Kurrer and Schulten proposed a moment method, expanding FPE for the $`N=\mathrm{}`$ model in a Taylor series around the center of distribution . Rodriguez and Tuckwell (RT) adopted a different moment method in which the original stochastic DEs are expanded in a series of fluctuations around means of state variables -. RT’s original moment method was first applied to Fitzhugh-Nagumo (FN) neuron model  and then FN neuron networks . RT’s moment method takes into account means, variances and covaricance of local variables, then the number of DEs for $`N`$-unit FN neuron networks is $`N_{eq}=N(2N+3)`$, which is, for example, 20300 for $`N=100`$. When we apply RT’s moment method to the $`N`$-unit AR network under consideration, we get $`N_{eq}=(1/2)N(N+3)`$, which is 5150 for $`N=100`$. This exponential increase in the number of DEs prevents us from calculations for a system with a realistic size.

Quite recently the present author has proposed an alternative moment method, which is hereafter referred to as a dynamical mean-field approximation (DMA) . In DMA we take into account means, variances and covariances of local and global variables, replacing the original $`2N`$-dimensional DEs of the $`N`$-unit FN model by eight-dimensional DEs independently of $`N`$. We have investigated the $`N`$ dependence of the spike timing precision, which has been shown to be improved by increasing $`N`$. The purpose of the present paper is to apply DMA to coupled AR networks to discuss their dynamics. Although in our previous paper , we included up to the forth-order moments, we will, in this paper, take into account up to infinite-order moments, which are expressed as a product of second-order moments with the use of the Gaussian assumption for the distribution of state variables. The number of equations for $`N`$-unit AR networks becomes $`N_{eq}=3`$ in DMA, which is much smaller than $`N_{eq}=(1/2)N(N+3)`$ in the conventional moment method.

The paper is organized as follows: In Sec. II, we develop a DMA theory for an ensemble of $`N`$ systems, obtaining equations of motions of moments of local and global variables. In Sec. III, the phase diagram showing various states in coupled AR models is discussed. Discussions are given in Sec. IV, where we compare our DMA with the conventional moment method, showing that the former may be alternatively derived from the latter with a proper reduction of the number of variables with the mean-field approximation. The final Sec. V is devoted to conclusions.

## II Dynamical mean-field approximation

### A Basic formulation

We assume an ensemble consisting of coupled $`N`$-unit phase models subject to white noises, whose dynamics of the phase $`\varphi _i`$ (mod $`2\pi `$) of the $`i`$th system is described by nonlinear DEs given by

$`{\displaystyle \frac{d\varphi _i(t)}{dt}}`$ $`=`$ $`F(\varphi _i)+{\displaystyle \frac{w}{N}}{\displaystyle \underset{j}{}}G(\varphi _j\varphi _i)+\xi _i(t),\text{(}i=1N\text{)}`$ (1)
where explicit forms of $`F(x)`$ and $`G(x)`$ will be given shortly \[Eq. (14)\], $`w`$ denotes the coupling, and $`\xi _i(t)`$ is the independent Gaussian white noise with $`<\xi _i(t)>=0`$ and $`<\xi _i(t)\xi _j(t^{})>=2D\delta _{ij}\delta (tt^{})`$, the bracket $`<>`$ denoting the average over stochastic random variables \[see Eq. (4)\].

We will express these nonlinear DEs by moments of local and global variables of the ensemble. The global variable is defined by

$`\mathrm{\Phi }(t)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\varphi _i(t),`$ (2)
and their averages by

$`\mu (t)`$ $`=`$ $`<\mathrm{\Phi }(t)>,`$ (3)
where the bracket $`<>`$ denotes

$`<G(\mathit{\varphi })>`$ $`=`$ $`{\displaystyle \mathrm{}𝑑\mathit{\varphi }G(\mathit{\varphi },t)p(\mathit{\varphi },t)},`$ (4)
$`p`$($`\mathit{\varphi }`$,$`t`$) being the probability distribution function (pdf) for the $`N`$-dimensional stochastic variable $`\mathit{\varphi }`$=$`(\varphi _1,\mathrm{}.,\varphi _N)`$.

We express DEs given by Eq. (1) in terms of the deviations from their average given by

$`\delta \varphi _i(t)`$ $`=`$ $`\varphi _i(t)\mu (t),`$ (5)
to get variances between local and global variables, given by (the argument $`t`$ is hereafter neglected)

$`\gamma `$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}<\delta \varphi _i^2>,`$ (6)
$`\rho `$ $`=`$ $`<\delta \mathrm{\Phi }^2>={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}<\delta \varphi _i\delta \varphi _j>,`$ (7)
where $`\delta \mathrm{\Phi }=\mathrm{\Phi }(t)\mu (t)`$. It is noted that $`\gamma `$ expresses the spatial average of fluctuations in local variables of $`\varphi _i`$ while $`\rho `$ denotes fluctuations in a global variable of $`\mathrm{\Phi }`$.

We assume that the noise intensity $`D`$ is weak. This allows us to expand the right-hand side of Eq. (11) around the average $`\mu `$, to get

$`{\displaystyle \frac{d\varphi _i}{dt}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(\mathrm{})}}{\mathrm{}!}}(\delta \varphi _i)^{\mathrm{}}+{\displaystyle \frac{w}{N}}{\displaystyle \underset{j}{}}{\displaystyle \underset{\mathrm{}=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(\mathrm{})}}{\mathrm{}!}}(\delta \varphi _j\delta \varphi _i)^{\mathrm{}}+\xi _i(t),`$ (8)
from which we get

$`{\displaystyle \frac{d\mu }{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}<{\displaystyle \frac{d\varphi _i}{dt}}>`$ (9)
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}{\displaystyle \underset{\mathrm{}=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(\mathrm{})}}{\mathrm{}!}}<(\delta \varphi _i)^{\mathrm{}}>+{\displaystyle \frac{w}{N^2}}{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}{\displaystyle \underset{\mathrm{}=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(\mathrm{})}}{\mathrm{}!}}<(\delta \varphi _j\delta \varphi _i)^{\mathrm{}}>,`$ (10)
where $`F^{(\mathrm{})}=F^{(\mathrm{})}(\mu )`$ and $`G^{(\mathrm{})}=G^{(\mathrm{})}(0)`$.

We furthermore assume that pdf of state variables takes Gaussian form. Numerical simulations have shown that for weak noises, the distribution of $`\varphi (t)`$ in a single AR system nearly obeys the Gaussian distribution, although for strong noises, the distribution of $`\varphi (t)`$ deviates from the Gaussian . Similar behavior of the Gaussian distribution of state variables has been reported in FN and Hodgkin-Huxley neuron models . When we adopt the Gaussian decoupling approximation, averages higher than the second-order moments in Eq. (9) may be expressed in terms of the second moments given by

$`<\delta \varphi _1,\mathrm{},\delta \varphi _{\mathrm{}}>`$ $`=`$ $`{\displaystyle \underset{allparings}{}}\mathrm{\Pi }_{km}<\delta \varphi _k\delta \varphi _m>,\text{for even }\mathrm{},`$ (11)
$`=`$ $`0,\text{for odd }\mathrm{},`$ (12)
where the summation is performed for all $`(\mathrm{}1)(\mathrm{}3)\mathrm{}.31`$ combinations. After some manipulations with the use of the approximations mentioned above, we get equations of motions for $`\mu `$, $`\gamma `$ and $`\rho `$ given by (for details see Appendix A)

$`{\displaystyle \frac{d\mu }{dt}}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n)}}{n!}}({\displaystyle \frac{\gamma }{2}})^n+w{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n)}}{n!}}(\gamma \rho )^n,`$ (13)
$`{\displaystyle \frac{d\gamma }{dt}}`$ $`=`$ $`2\gamma {\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n+1)}}{n!}}({\displaystyle \frac{\gamma }{2}})^n+2w(\gamma \rho ){\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n+1)}}{n!}}(\gamma \rho )^n+2D,`$ (14)
$`{\displaystyle \frac{d\rho }{dt}}`$ $`=`$ $`2\rho {\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n+1)}}{n!}}({\displaystyle \frac{\gamma }{2}})^n+{\displaystyle \frac{2D}{N}}.`$ (15)
The coupled AR network is given by

$`{\displaystyle \frac{d\varphi _i(t)}{dt}}`$ $`=`$ $`casin(x)+{\displaystyle \frac{w}{N}}{\displaystyle \underset{j}{}}sin(\varphi _j\varphi _i)+\xi _i(t),\text{(}i=1N\text{)}`$ (16)
with $`F(x)=casin(x)`$ and $`G(x)=sin(x)`$ in Eq. (1). In Eq. (14) $`c(>0)`$ stands for the intrinsic frequency and $`a`$ the intensity of the pinning force introduced such that for $`c<a`$, the system mimics the stochastic limit cycle or excitable elements. The model with $`c=0`$ stands for the equilibrium planar model. The case of $`a=0`$ corresponds to a usual phase model . For $`w=0`$ and $`\xi _i=0`$, the AR system locates at the stationary point given by $`\varphi _i=\varphi ^{}=arcsin(c/a)`$. When noises are introduced, the system shows the intrigue behavior. Substituting $`F(x)`$ and $`G(x)`$ to Eqs. (11)-(13), we obtain DEs for $`\mu `$, $`\gamma `$ and $`\rho `$ given by

$`{\displaystyle \frac{d\mu }{dt}}`$ $`=`$ $`casin(\mu )exp({\displaystyle \frac{\gamma }{2}}),`$ (17)
$`{\displaystyle \frac{d\gamma }{dt}}`$ $`=`$ $`2a\gamma cos(\mu )exp({\displaystyle \frac{\gamma }{2}})2w(\gamma \rho )exp[(\gamma \rho )]+2D,`$ (18)
$`{\displaystyle \frac{d\rho }{dt}}`$ $`=`$ $`2a\rho cos(\mu )exp({\displaystyle \frac{\gamma }{2}})+{\displaystyle \frac{2D}{N}},`$ (19)
The original $`N`$-dimensional stochastic DEs given by Eq. (1) are transformed to three-dimensional deterministic DEs, which show much variety depending on model parameters such as $`a`$, $`c`$, $`w`$, $`D`$ and $`N`$.

We note that the noise contribution is $`2D`$ in Eq. (12) while that is $`2D/N`$ in Eq. (13). It is easy to see that

$`\rho `$ $`=`$ $`\gamma /N,\text{(for }w/D0\text{)}`$ (20)
$`=`$ $`\gamma .\text{(for }D/w0\text{)}`$ (21)
Equation (18) agrees with the central-limit theorem. In the limit of $`N=\mathrm{}`$, we get $`\rho =0`$. On the contrary, in the limit of $`N=1`$, we have $`\rho =\gamma `$.

### B Various quantities

Distribution of local variables

Adopting the mean-field approximation, we get $`<\varphi _i>(1/N)_k<\varphi _k>=\mu `$ and $`<\delta \varphi _i^2>(1/N)_k<\delta \varphi _k^2>=\gamma `$. Then the distribution for the variable $`\varphi _i`$ is given by

$`P(\varphi _i)`$ $``$ $`({\displaystyle \frac{1}{\sqrt{\gamma }}})\overline{\varphi }({\displaystyle \frac{\varphi _i\mu }{\sqrt{\gamma }}}),`$ (22)
where $`\overline{\varphi }(x)`$ is the normal distribution function given by

$$\overline{\varphi }(x)=\frac{1}{\sqrt{2\pi }}\mathrm{exp}(\frac{x^2}{2}).$$
(23)
The probability given by Eq. (20) depends on the time because of the time dependence of $`\mu (t)`$ and $`\gamma (t)`$.

Distribution of global variables

Mean and variance of global variables $`\mathrm{\Phi }`$ are given by $`<\mathrm{\Phi }>=\mu `$ and $`<\delta \mathrm{\Phi }^2>=\rho `$, respectively. We get the distribution for the global variable $`\mathrm{\Phi }`$ given by

$`P(\mathrm{\Phi })`$ $``$ $`({\displaystyle \frac{1}{\sqrt{\rho }}})\overline{\varphi }({\displaystyle \frac{\mathrm{\Phi }\mu }{\sqrt{\rho }}}).`$ (24)
Averaged frequency

The averaged frequency $`\nu `$ is defined by

$$\nu =[\frac{1}{N(N_{fi}1)}\underset{i}{}\underset{k}{}T_{oik}]^1,$$
(25)
with

$$T_{oik}=t_{ik+1}t_{ik},$$
(26)
$$t_{ik}=\{t\varphi _i(t)=\theta ;\dot{\varphi _i}>0;tt_{ik1}+\tau _r\},$$
(27)
where $`N_{fi}`$ stands for the number of firings of a given rotator $`i`$, $`\dot{\varphi _i}=d\varphi _i/dt`$, $`T_{oik}`$ the interspike interval (ISI) of output signals, $`t_{ik}`$ the $`k`$th firing time, and $`\theta (=2\pi )`$ and $`\tau _r(=5)`$ are the threshold level and the refractory period, respectively. When there is no firings, we set $`\nu =0`$.

Order parameters

The order parameter $`\zeta `$ and its fluctuation $`\delta \zeta `$ are defined by

$`\zeta `$ $`=`$ $`\overline{z(t)},`$ (28)
$`\delta \zeta `$ $`=`$ $`\sqrt{\overline{z(t)^2}\zeta ^2},`$ (29)
with

$`z(t)`$ $`=`$ $`{\displaystyle \underset{i}{}}exp[i\varphi _i(t)]=<exp[i\varphi _i(t)]>,`$ (30)
where the overline denotes the temporal avarage. By expanding $`z(t)`$ in a series of $`\delta \varphi _i`$ around $`\mu `$ and adopting the Gaussian decoupling approximation given by Eq. (10), we get

$$z(t)=exp[i\mu (t)]\underset{n}{}\frac{1}{n!}(\frac{\gamma }{2})^n=exp[i\mu (t)\frac{\gamma (t)}{2}]$$
(31)
Synchronization ratio

The synchronization ratio $`\sigma `$ is defined by

$$\sigma =\overline{s(t)},$$
(32)
with

$$s(t)=\frac{(\rho /\gamma 1/N)}{(11/N)}.$$
(33)
For the completely synchronous (asynchronous) state, both $`\zeta `$ and $`\sigma `$ become 1 (0). It is noted, however, that while $`\zeta `$ depends on $`\gamma `$, $`\sigma `$ is a function of $`(\rho /\gamma 1/N)=N^1_{ij}<\delta \varphi _i\delta \varphi _j>/_i<\delta \varphi _i^2>`$; the ratio of the inter-AR correlation to the intra-AR correlation.

Before discussing calculated results with the use of DMA, it is worth to mention the calculation of Kurrer and Schulten . They intended to solve the FPE for the $`N=\mathrm{}`$ AR model given by

$$\frac{}{t}n(\varphi ,t)=\frac{}{\varphi }[F(\varphi ,t)w𝑑\varphi ^{}sin(\varphi \varphi ^{})n(\varphi ^{},t)]n(\varphi ,t)+D\frac{^2}{t^2}n(\varphi ,t),$$
(34)
where the density probability $`n(\varphi ,t)`$ is defined by

$`n(\varphi ,t)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\delta (\varphi \varphi _i(t)),`$ (35)
with the periodic condition: $`n(\varphi +2\pi ,t)=n(\varphi ,t)`$ and the normalization condition: $`𝑑\varphi n(\varphi ,t)=1`$. Kurrer and Schulten  expanded $`F(\varphi ,t)`$ in a Taylor series around the center of distribution, assuming the Gaussian form for $`n(\varphi ,t)`$ given by

$`n(\varphi ,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi v(t)}}}exp\left({\displaystyle \frac{[\varphi u(t)]^2}{2v(t)}}\right),`$ (36)
where the mean $`u(t)`$ and variance $`v(t)`$ obey DEs given by

$`{\displaystyle \frac{du}{dt}}`$ $`=`$ $`casin(u)exp({\displaystyle \frac{v}{2}}),`$ (37)
$`{\displaystyle \frac{dv}{dt}}`$ $`=`$ $`2avcos(u)exp({\displaystyle \frac{v}{2}})2wvexp(v)+2D.`$ (38)
Equations (35) and (36) resemble our Eqs. (15)-(17) if we read $`u\mu `$ and $`v\gamma `$. Actually, Eqs. (35) and (36) are equvalent to Eqs. (15) and (16) in the case of $`\rho =0`$, which is realized in the limit of $`N=\mathrm{}`$.

## III Calculated results

DMA equations given by Eqs. (15)-(17) have been solved by the forth-order Runge-Kutta method with a time step of 0.01, the initial conditions being $`\mu (0)=\gamma (0)=\rho (0)=0`$. Calculations have been performed for $`0t1000`$ (100 000 steps) and results for $`t<100`$ are discarded. Simulations of directly solving Eq. (1) have been made by the forth-order Runge-Kutta method with a time step of 0.01, the initial conditions being $`\varphi _i(0)=0`$ ($`i=1`$ to $`N`$). The number of trials for a given set of parameters in our simulations is hundred otherwise noticed. We have solved also FPE given by Eq. (32), which is valid for the $`N=\mathrm{}`$ AR model. We first Fourier transform FPE with the first 30 modes after Ref. . A set of 61 ordinary DEs has been solved by the Runge-Kutta method.

### A Phase diagram for various types of solutions

By solving Eqs. (15)-(17), we get the stationary state and the non-stationary state: in the former state the variables are time independent while in the latter state they are time dependent. The equilibrium values of $`\mu `$, $`\gamma `$ and $`\rho `$ in the stationary state are given by (we set $`c=1`$ hereafter)

$`\mu `$ $`=`$ $`arcsin[({\displaystyle \frac{1}{a}})exp({\displaystyle \frac{\gamma }{2}})],`$ (39)
$`\gamma `$ $`=`$ $`{\displaystyle \frac{D+w\rho exp[(\gamma \rho )]}{\sqrt{a^2exp(\gamma )1}+wexp[(\gamma \rho )]}},`$ (40)
$`\rho `$ $`=`$ $`{\displaystyle \frac{D/N}{\sqrt{a^2exp(\gamma )1}}}.`$ (41)
The stationary state where Eqs. (37)-(39) are satisfied, is hereafter referred to as the S state. Kurrer and Schulten  pointed out that the non-stationary state may be classified into the time periodic (P) state and the random, disordered (R) state. DMA also yields three types of the S, P and R states characterized by the quantities of $`\zeta `$, $`\delta \zeta `$, $`\nu `$ and $`\sigma `$ introduced in Sec. IIB, result being summarized in Table 1. In the S state, $`\delta \zeta `$ and $`\nu `$ are vanishing while $`\zeta (1)`$ and $`\sigma `$ are finite. In the P state, all quantities are finite. In contrast, in the R state, all quantities except $`\nu `$ vanish.

Boundaries between these three states depend on $`a`$, $`w`$, $`D`$ and $`N`$. Figure 1 expresses the $`Da`$ phase diagram showing the boundaries between these states in coupled AR models calculated with the use of DMA for $`N=10`$, 100 and $`N=\mathrm{}`$. The gradient of the boundary between the stationary (S) state and non-stationary (P+R) states is decreased as increasing the value of $`w`$ and/or of $`N`$. The difference between boundaries for $`N=10`$ and $`N=\mathrm{}`$ with $`w=0.1`$ is very small: the effect of $`N`$ becomes more significant for a stronger coupling. The critical $`a`$ value, $`a_c`$, above which the S state exists, is given by

$$a_c1[c_1c_2w(1\frac{1}{N})]D,$$
(42)
where $`c_1=2.25`$ and $`c_2=1.75`$. In contrast, the critical value of $`a_d`$ for the boundary between the P and R states for $`w=1.0`$ is given by

$$a_d1(d_1+\frac{d_2}{N})[D(d_3+\frac{d_4}{N})],$$
(43)
where $`d_1=5.36`$, $`d_2=257`$, $`d_3=0.265`$ and $`d_4=0.9`$.

The behavior of the solutions of DMA in the S, P and R states when $`D`$ and/or $`a`$ values are changed, is shown in Figs. 2-5. We will first mention the calculations of DMA in the three states. Figure 2(a) and 2(b) express the $`D`$ dependence of $`\zeta `$, $`\delta \zeta `$, $`\nu `$ and $`\sigma `$ for a representative set of parameters of $`a=1.05`$, $`w=1.0`$ and $`N=100`$, showing that the network is in the S state for $`D0.082`$, in the P state for $`0.082<D0.273`$, and in the R state for $`D>0.273`$. In contrast, Fig. 3(a) and 3(b) express $`\zeta `$, $`\delta \zeta `$, $`\nu `$ and $`\sigma `$ as a function of $`a`$ for a set of parameters of $`D=0.1`$, $`w=1.0`$ and $`N=100`$, for which the networks is in the S state for $`a1.06`$ and in the P state for $`a<1.06`$.

Equations (37)-(39) for small $`D`$ and $`w`$ in the S state yield

$`\mu `$ $`=`$ $`arcsin({\displaystyle \frac{1}{a}})+d_1D+..,`$ (44)
$`\gamma `$ $`=`$ $`{\displaystyle \frac{D}{\sqrt{a^21}}}[1+d_2Dd_3(1{\displaystyle \frac{1}{N}})w+d_4w^2+..],`$ (45)
$`\rho `$ $`=`$ $`{\displaystyle \frac{(D/N)}{\sqrt{a^21}}}[1+d_2D+..,],`$ (46)
leading to

$`\zeta `$ $`=`$ $`1{\displaystyle \frac{\gamma }{2}},`$ (47)
$`\sigma `$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \frac{d_3w}{(1+d_2D)}},`$ (48)
where $`d_1=1/2(a^21)`$, $`d_2=a^2/2(a^21)^{3/2}`$, $`d_3=1/\sqrt{a^21}`$, and $`d_4(>\mathrm{\hspace{0.33em}0})`$ is a complex function of $`D`$, $`w`$ and $`N`$.

Solid curves in Figs. 4(a) and 4(b) show distributions of local \[$`P(\varphi _i(t))`$\] and global variables \[$`P(\mathrm{\Phi }(t))`$\], respectively, in DMA for $`a=1.05`$, $`w=1.0`$, $`D=0.05`$. They are obtained by Eqs. (20) and (22) with $`\mu =1.339`$, $`\gamma =0.04354`$ and $`\rho =0.00212`$.

Figures 2(a) and 2(b) show that when the noise intensity is increased and crosses the value of 0.082, the AR network begins correlated firings. This implies the appearance of the P state, where $`\delta \zeta `$, $`\nu `$ and $`\sigma `$ are continuously changed. Solid curves in Fig. 4(c) ad 4(d) express $`P(\varphi _i(t))`$ and $`P(\mathrm{\Phi }(t))`$ for $`D=0.10`$, respectively, which are given by Eqs. (20) and (22) with $`\mu =1.497`$, $`\gamma =0.11022`$ and $`\rho =0.009443`$. The time evolution of the probability of $`P(\varphi (t))`$ calculated in DMA for $`D=0.10`$ in the P state is shown in Fig. 5(a), which is oscillating with the period of about 40. It is noted that not only the position of $`P(\varphi (t))`$ but also its width change as a function of $`t`$. For example, we get $`\mu =6.151`$ and $`\gamma =1.711`$ at $`t=120`$ while $`\mu =1.497`$ and $`\gamma =0.11022`$ at $`t=100`$.

When the $`D`$ value is rmore increased, the AR network fires abundantly and irregularly, which suggests the appearance of the R state. The solution of Eqs. (37)-(39) in the R state for a large $`t`$ is given by

$`\mu `$ $``$ $`ct,`$ (49)
$`\gamma `$ $``$ $`2Dt,`$ (50)
$`\rho `$ $``$ $`({\displaystyle \frac{2D}{N}})t,`$ (51)
which lead to vanishing $`\zeta `$, $`\delta \zeta `$ and $`\sigma `$ except $`\nu `$. Figure 2(a) and 2(b) show that $`\zeta `$ and $`\delta \zeta `$ suddenly vanish at $`D=0.273`$ with no hysteresis. Solid curves in Fig. 4(e) ad 4(f) express $`P(\varphi _i(t))`$ and $`P(\mathrm{\Phi }(t))`$ for $`D=0.30`$, respectively, which are given by Eqs. (20) and (22) with $`\mu =1.339`$, $`\gamma =19.9`$ and $`\rho =0.199`$.

In the following, results of DMA will be compared with those of simulations and FPE. Dashed curves in Figs. 2(a) and 2(b) show the results of simulations for the $`D`$ dependence of $`\zeta `$, $`\delta \zeta `$, $`\nu `$ and $`\sigma `$. The agreement of $`\zeta `$ in DMA with that in simulations is good for S and P states. However, it is not good in the R state, where $`\zeta `$ vanishes in DMA but not in simulations. This is expected to be due to deviations of the state-variable distributions from the Gaussian form. When $`D`$ is more increased in the R state, our simulations yields a gradual decrease in $`\zeta `$, which is 0.454, 0.276, 0.188, 0.139 and 0.110 for $`D=1.0`$, 2.0, 3.0, 4.0 and 5.0, respectively, with $`a=1.05`$, $`w=1.0`$ and $`N=100`$. We note in Fig. 2(a) that $`\delta \zeta `$ of simulations is about ten times smaller than that of DMA. This is clearly seen in Fig. 6(a) where we plot the time evolution of $`z(t)`$ obtained by DMA and simulations for $`a=1.05`$, $`w=1.0`$, $`D=0.10`$ and $`N=100`$. The former has larger temporal fluctuations than the latter although both yield similar averaged values of $`\zeta =\overline{z(t)}`$. In contrast, Fig. 6(b) shows the time dependence of $`s(t)`$ calculated by DMA and simulations for $`a=1.05`$, $`w=1.0`$, $`D=0.10`$ and $`N=100`$. Again our DMA yields larger fluctuations in $`s(t)`$ than simulations although both methods lead to similar averaged values of $`\sigma =\overline{s(t)}`$. When comparing K(a) with K(b), we notice that the time dependence of $`z(t)`$ is not the same as that of $`s(t)`$. This is because $`z(t)`$ is a function of $`\gamma `$ \[Eq. (29)\] while $`s(t)`$ is a function of the ratio of $`\rho /\gamma `$ \[Eq. (31)\]. Dashed curves in Figs. 3(a) and 3(b) show the $`a`$ dependence of $`\zeta `$, $`\delta \zeta `$, $`\nu `$ and $`\sigma `$ obtained by simulations. $`P(\varphi )`$ and $`P(\mathrm{\Phi })`$ obtained by simulations are plotted by dashed curves in Fig. 4(a)-4(f). Dotted curves in Fig. 4(a), 4(c) and 4(e) denote $`n(\varphi )`$ obtained by FPE for the $`N=\mathrm{}`$ AR model. Figure 5(b) express the time evolution of $`n(\varphi ,t)`$ in the P state obtained by FPE. From a comparison between the results of DMA and simulations (and FPE) mentioned above, we note that DMA is good for the S state, in fairly good for the P state in the qualitative sense, but not good for the R state.

For a comparison, we show by the dotted curve in Fig. 1, the boundary obtained by Shimokawa and Kuramoto (SK) with the use of the FPE for $`w=1.0`$ and $`N=\mathrm{}`$ . The ordered P state where $`\delta \zeta 0`$ and $`\nu 0`$ is reported to exist in the triangle region enclosed by the dotted curve and the horizontal axis. The P state obtained by SK is nearly in agreement with our P state. In SK, states besides the P state are regarded as the stationary state where $`n(\varphi ,t)/t=0`$. On the contrary, Kurrer and Schulten (KS) distinguished the R state from the P state, both of which are non-stationary ($`\nu 0`$) . The results of SK and KS are for the $`N=\mathrm{}`$ AR model. Figures 7(a) and 7(b) show the $`D`$ dependence of $`\zeta `$, $`\delta \zeta `$ and $`\nu `$ for a set of parameters of $`a=1.05`$, $`w=1.0`$ and $`N=10^4`$, which are the same as in Figs. 5(a) and 5(b) except $`N`$. In order to simulate the $`N=\mathrm{}`$ limit, the $`N`$ value in Figs. 7(a) and 7(b) is chosen to be very large but finite because $`s(t)`$ given by Eq. (31) is not properly defined in this limit. Results for $`N=\mathrm{}`$ in Figs. 1 and 7 should be compared with those obtained by SK and KS. As was pointed in Sec. IIB, DEs of DMA given by Eqs. (15)-(17) in the limit of $`N=\mathrm{}`$ agree with those of KS given by Eqs. (35)-(36). Nethertheless, our $`Da`$ phase diagram for $`N=\mathrm{}`$ in Fig. 1 does not agree with that of KS. For example, KS obtained the critical values given by

$`a_c1`$ $`=`$ $`{\displaystyle \frac{D}{2w}},\text{(between S and P + R)}`$ (52)
$`{\displaystyle \frac{D_d}{w}}`$ $`=`$ $`0.736,\text{(between P and R)}`$ (53)
which do not agree with our expressions given by Eqs. (40) and (41) for $`N=\mathrm{}`$.

### B Cluster-size ($`N`$) dependence in the S state

Since one of the advantages of DMA is that we can discuss the finite-$`N`$ property of coupled AR networks, we have made more detailed calculations of the $`N`$ dependence of the quantities in the S state. Figures 8(a) and 8(b) show the log-log plot of the $`N`$ dependence of $`\gamma `$, $`\rho `$ and $`\sigma `$, results for $`N=10`$ and 20 being for 500 trials. Solid curves in Fig. 8(a) express the result of DMA for a set of parameters of $`a=1.05`$, $`w=1.0`$ and $`D=0.05`$, whereas circles, squares and triangles denote those of simulations. For this set of parameters, the P state is realized for $`N9`$. We note that as $`N`$ is decreased from above and approaches to the S-P boundary, fluctuations of $`\gamma `$ and $`\rho `$ are increased (and $`\sigma `$ is also increased). Similar behavior is observed for a different set of parameters. Figure 8(b) shows the results for $`a=1.20`$, $`w=1.0`$ and $`D=0.10`$. With this set of parameters, we get the S state for $`N2`$ and the P state for $`N=1`$. Figures 8(a) and 8(b) show that as increasing $`N`$, $`\rho `$ is much decreased but $`\gamma `$ shows only a weak $`N`$ dependence. We note that $`\sigma `$ is decreased as increaing $`N`$ whereas $`\zeta =exp(\gamma /2)`$ shows little $`N`$ dependence. We should stress that although fluctuations of global variables is inversely decreased as $`\rho 1/N`$ consistent with the central-limit theorem, those of local variables remain finite even for $`N\mathrm{}`$.

## IV Discussions

We have proposed DMA theory for stochastic, nonlinear networks like coupled AR models, taking into account means, variances and covariances of local and global variables. It is worth to compare DMA with the conventional moment method in which means, variances and covariances of local variables are given by

$`m_i`$ $`=`$ $`<\varphi _i>,`$ (54)
$`C_{ij}`$ $`=`$ $`<\mathrm{\Delta }\varphi _i\mathrm{\Delta }\varphi _j>,`$ (55)
with $`\mathrm{\Delta }\varphi _i=\varphi _im_i`$. By using the Gaussian decoupling approximation \[Eq. (10)\], we get (for details see Appendix B)

$`{\displaystyle \frac{dm_i}{dt}}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F_i^{(2n)}}{2^nn!}}C_{ii}^n+{\displaystyle \frac{w}{N}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n)}}{2^nn!}}[C_{kk}+C_{ii}2C_{ik}]^n,`$ (56)
$`{\displaystyle \frac{dC_{ij}}{dt}}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F_i^{(2n+1)}}{2^nn!}}(C_{jj}^n+C_{ii}^n)C_{ij}+{\displaystyle \frac{w}{N}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n+1)}}{2^nn!}}`$ (58)
$`\times [(C_{kk}+C_{ii}2C_{ik})^n(C_{jk}C_{ij})+(C_{kk}+C_{jj}2C_{jk})^n(C_{ik}C_{ij})]+2D_i\delta _{ij},`$
where $`F_i^{(\mathrm{})}=F^{(\mathrm{})}(m_i)`$ and $`G^{(\mathrm{})}=G^{(\mathrm{})}(0)`$.

For the AR model, Eqs. (54) and (55) become

$`{\displaystyle \frac{dm_i}{dt}}`$ $`=`$ $`casin(m_i)exp({\displaystyle \frac{1}{2}}C_{ii}),`$ (59)
$`{\displaystyle \frac{dC_{ij}}{dt}}`$ $`=`$ $`acos(m_i){\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^n}{2^nn!}}(C_{ii}^n+C_{jj}^n)C_{ij}+2D\delta _{ij}`$ (62)
$`+{\displaystyle \frac{w}{N}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^n}{2^nn!}}[(C_{ii}+C_{kk}2C_{ik})^n(C_{jk}C_{ij})`$
$`+(C_{jj}+C_{kk}2C_{jk})^n(C_{ik}C_{ij})],`$
For variances ($`i=j`$), Eq. (57) becomes

$`{\displaystyle \frac{dC_{ii}}{dt}}`$ $`=`$ $`2acos(m_i)C_{ii}exp({\displaystyle \frac{1}{2}}C_{ii})+2D`$ (64)
$`+{\displaystyle \frac{2w}{N}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^n}{2^nn!}}[(C_{ii}+C_{kk}2C_{ik})^n(C_{ik}C_{ii})].`$
Taking into the symmetry relations: $`C_{ij}=C_{ji}`$, we get the number of DEs in the moment method to be $`N_{eq}=N(N+3)/2`$, which is 65, 5150 and 501 500 for $`N=`$10, 100 and 1000, respectively, while $`N_{eq}=3`$ in our DMA.

It will be shown that we can derive DMA from the moment method by reducing the number of DEs, adopting the mean-field approximation:

$`m_i`$ $``$ $`\mu ,`$ (65)
$`C_{ii}^n`$ $``$ $`\gamma ^{n1}C_{ii},`$ (66)
$`(C_{kk}+C_{ii}2C_{ik})^n`$ $``$ $`2^{n1}(\gamma \rho )^{n1}(C_{kk}+C_{ii}2C_{ik}),\text{(}ik\text{)}`$ (67)
with the relations given by

$`\mu ={\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}m_i,`$ (68)
$`\gamma ={\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}C_{ii},`$ (69)
$`\rho ={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}C_{ij}.`$ (70)
DEs for $`\mu `$, $`\gamma `$ and $`\rho `$ are given by Eqs.(11)-(13) for the general phase model or by Eqs. (15)-(17) for the AR model.

It is possible to regard DMA as the single-site self-consistent theory. Let us assume a configuration in which a single nonlinear system $`i`$ is embedded in an effective medium whose effect is realized by a given system $`i`$ as its effective external input through the coupling $`w`$. Assuming $`m_i=\mu `$, we replace quantities in coupling terms of Eqs. (57) and (58) by effective quantities of $`\mu `$, $`\gamma `$ and $`\rho `$, to get

$`{\displaystyle \frac{dm_i}{dt}}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n)}}{2^nn!}}C_{ii}^n+{\displaystyle \frac{w}{N}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n)}}{n!}}(\gamma \rho )^n,`$ (71)
$`{\displaystyle \frac{dC_{ij}}{dt}}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n+1)}}{2^nn!}}(C_{jj}^n+C_{ii}^n)^nC_{ij}+{\displaystyle \frac{w}{N}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n+1)}}{n!}}(\gamma \rho )^{n+1}.`$ (72)
We should note that $`m_i`$ and $`C_{ij}`$ determined by Eqs. (65) and (66) are functions of $`\mu `$, $`\gamma `$ and $`\rho `$. In order to self-consistently determine them, we impose the self-consistent conditions given by

$`\mu `$ $`=`$ $`m_i,`$ (73)
$`\gamma `$ $`=`$ $`C_{ii},`$ (74)
$`\rho `$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{j}{}}C_{ij}.`$ (75)
Note that Eqs. (67)-(69) are assumed to hold independently of $`i`$ and that $`m_i`$ and $`C_{ij}`$ in their right-hand sides are functions of $`\gamma `$ and $`\rho `$. The condition given by Eqs. (65)-(69) with the mean-field approximation given by Eq. (59)-(61) yields DEs for $`\gamma `$ and $`\rho `$ which are again given by Eqs. (11)-(13). The self-consistent condition given by Eq. (67)-(69), which assumes that the quantities averaged at a given site are the same as those of the effective medium, is common in mean–field theories such as the Weiss theory for magnetism  and the coherent-potential approximation for random alloys .

By using DMA, we have investigated the response of the excitable, coupled AR networks to an applied periodic spike, by adding to the right-hand side of Eq. (15), the input term given by

$`I_{in}(t)`$ $`=`$ $`g,\text{for }mT_pt<mT_p+T_w\text{ (m: integer)}`$ (76)
$`=`$ $`0,\text{otherwise}`$ (77)
where $`g`$ denotes the magnitude, and $`T_p`$ (=50) and $`T_w`$ (=5) stand for the period and the duration of spikes, respectively. We get the critical value of $`g_c=0.159`$ below which there are no firings for $`D=0`$. Figure 9(a) shows the distribution of ISI, $`T_o`$, of output signals defined by Eqs. (24) and (25) as a function of $`D`$ in the absence of input spikes ($`g=0`$) for $`a=1.05`$, $`w=1.0`$ and $`N=100`$. Firings begin at $`D=0.082`$, above which the system is in the P state as discussed in Sec. IIIA. Around the P-R transition at $`D=0.273`$, ISIs have a small distribution. When the input spike is applied, distributions of ISIs are significantly changed. Figure 9(b) shows the distribution when the subthreshold input with $`g=0.1`$ $`(<g_c)`$ is applied. Firings occur at $`D0.04`$ with a help of noises. A flat segment at $`0.04<D<0.08`$ corresponds to a periodic solution locked to input spikes while the others show the complex behavior. In contrast, Fig. 9(c) shows the distribution of ISIs for the suprathreshold input with $`g=0.2`$. At $`0.05<D<0.08`$ in the S state, a new branch with $`35<T_o<43`$ appears beside the branch with $`T_o=50`$ locked to inputs. The distribution of ISIs in the presence of input spikes has much variety than that in the absence of noises, in particular in the P state, where the bifurcation is realized as the noise intensity is changed.

It is possible to discuss the firing-time accuracy of excitable AR models for an external input with the use of DMA . The $`k`$th firng time of a given rotator $`i`$ is defined as the time when $`\varphi _i(t)`$ crosses the threshold $`\theta `$ from below \[Eq. (25)\]:

$$t_{ik}=\{t\varphi _i(t)=\theta ;\dot{\varphi _i}>0;tt_{ik1}+\tau _r\}.$$
(78)
By using the discussion presented in Sec. IIB, the probability $`W_{\mathrm{}}`$ when $`\varphi _i(t)`$ at $`t`$ is above the threshold $`\theta `$ is given by

$$W_{\mathrm{}}(t)=1\psi (\frac{\theta \mu }{\sqrt{\gamma }}),$$
(79)
where $`\psi (y)`$ is the error function given by integrating the normal distribution function $`\overline{\varphi }(x)`$ from $`\mathrm{}`$ to $`y`$ \[Eq. (21)\]. The fraction of a given rotator $`i`$ emitting output at $`t`$ is given by

$$Z_{\mathrm{}}(t)=\frac{dW_{\mathrm{}}}{dt}\mathrm{\Theta }(\dot{\mu })=\overline{\varphi }(\frac{\theta \mu }{\sqrt{\gamma }})\frac{d}{dt}(\frac{\mu }{\sqrt{\gamma }})\mathrm{\Theta }(\dot{\mu }),$$
(80)
where $`\mathrm{\Theta }(x)=1`$ for $`x0`$ and 0 otherwise, and $`\dot{\mu }=d\mu (t)/dt`$. When we expand $`\mu (t)`$ in Eq. (73) around $`t_o^{}`$ where $`\mu (t_o^{})=\theta `$, we get

$$Z_{\mathrm{}}(t)\overline{\varphi }(\frac{tt_o^{}}{\delta t_o\mathrm{}})\frac{d}{dt}(\frac{\mu }{\sqrt{\gamma }})\mathrm{\Theta }(\dot{\mu }),$$
(81)
with

$$\delta t_o\mathrm{}=\frac{\sqrt{\gamma }}{\dot{\mu }}.$$
(82)
We note that $`Z_{\mathrm{}}`$ provides the distribution of firing times, showing that most of firings locate in the range given by

$$t_o\mathrm{}[t_o^{}\delta t_o\mathrm{},t_o^{}+\delta t_o\mathrm{}].$$
(83)
In the limit of vanishing $`D`$, Eq. (74) reduces to

$$Z_{\mathrm{}}(t)=\delta (tt_o^{}).$$
(84)
Similarly, we define the $`k`$th firing time relevant to the global variable $`\mathrm{\Phi }(t)`$ as

$$t_{gk}=\{t\mathrm{\Phi }_i(t)=\theta ;\dot{\mathrm{\Phi }_i}>0;tt_{gk1}+\tau _r\}.$$
(85)
The distribution of firing times $`t_g`$ is given by

$$Z_g(t)\overline{\varphi }(\frac{tt_o^{}}{\delta t_{og}})\frac{d}{dt}(\frac{\mu }{\sqrt{\rho }})\mathrm{\Theta }(\dot{\mu }),$$
(86)
with

$$\delta t_{og}=\frac{\sqrt{\rho }}{\dot{\mu }}.$$
(87)
Equation (79) shows that most of $`t_{og}`$ locate in the range given by

$$t_{og}[t_o^{}\delta t_{og},t_o^{}+\delta t_{og}].$$
(88)
From Eqs. (75) and (80), we get

$$\frac{t_{og}}{t_o\mathrm{}}=\sqrt{\frac{\rho }{\gamma }}\frac{1}{\sqrt{N}}.\text{(as }w/D0\text{)}$$
(89)
This implies that the firing-time accuracy is improved as the ensemble size is increased even when there no couplings among ARs. This is consistent with results reported prevously -.

## V Conclusions

We have developed DMA, which has been shown to be derived in various ways: equations of motions for means, variances and covariances of local and global variables (Sec. IIA), a reduction in the number of moments in the moment method, and a single-site self-consistent approximation to the moment method (Sec. IV). Our DMA theory, which assumes weak noises and the Gaussian distribution of state variables, goes beyond the weak coupling because no constraints are imposed on the coupling strength. The advantage of DMA is that it can be applied to large but finite-$`N`$ nonlinear systems subject not only to white noises but also to color noises. This is in contrast with FPE, which is applicable to a single or infinite system subject to white noises. The limitation of DMA is the weak noise, for which calculated results based on DMA are in fairly good agreement with those obtained by direct simulations. When the noise intensity becomes stronger, the state-variable distribution more deviates from the Gaussian form, and the agreement of results of DMA with those of simulations becomes worse. Nevertheless, our DMA is expected to be meaningful for qualitative or semi-quantitative discussion on the properties of coupled nonlinear systems. It is possible to regard DEs given by Eqs. (15)-(17) as the mean-filed AR model which may show interesting behavior for applied input signals and noises. When we consider an ensemble of $`N`$-unit systems, each of which is described by a $`M`$-variable nonlinear DE, the number of the deterministic DEs is $`N_{eq}=M+M(M+1)=M(M+2)`$ independently of $`N`$ in DMA while it is $`N_{eq}=NM+(1/2)NM(NM+1)=(1/2)NM(NM+3)`$ in the conventional moment. In the case of $`M=2`$ (as in FN model), for example, DMA leads to $`N_{eq}=8`$ while the moment method yields $`N_{eq}=N(2N+3)`$, which is 2310, 20 300 and 2 003 000 for $`N=10`$, 100 and 1000, respectively. These figures clearly show the advantage and feasibility of our DMA theory.

To summarized, the property of excitable AR networks has been discussed with the use of DMA. The obtained results are summarized as follows. (1) Depending on model parameters of $`a`$, $`w`$, $`D`$ and $`N`$, AR networks display three types of dynamics (Fig. 1): S, P and R states are characterized by the quantities of $`\zeta `$, $`\delta \zeta `$, $`\nu `$ and $`\sigma `$, as summarized in Table 1. (2) The S-P transition is of the (continuous) second-order one while P-R and S-R transitions are of the (discontinuous) first-order one with no hysteresis. (3) There are no enhancements in order-parameter fluctuations of $`\delta \zeta `$ at the transitions. (4) Fluctuations in local variables ($`\gamma `$) remain finite even for $`N=\mathrm{}`$ whereas those ($`\rho `$) in global variables varies as $`\rho 1/N`$, which is consistent with the central-limit theorem.

## Acknowledgements

This work is partly supported by a Grant-in-Aid for Scientific Research from the Japanese Ministry of Education, Culture, Sports, Science and Technology.

## A Derivation of Eqs. (11)-(13)

When we adopt the Gaussian decoupling approximation given by Eqs. (10), Eq. (9) becomes

$`{\displaystyle \frac{d\mu }{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n)}}{(2n)!}}B_{2n}<\delta \varphi _i^2>^n`$ (A2)
$`+{\displaystyle \frac{w}{N^2}}{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n)}}{(2n)!}}B_{2n}<(\delta \varphi _j\delta \varphi _i)^2>^n,`$
where $`B_{2n}=(2n1)(2n3)\mathrm{}.31`$. In deriving Eq. (A1), we treat $`(\delta \varphi _j\delta \varphi _j)`$ as a new variable with the Gaussian distribution. Adopting the mean-field approximation given by

$`<\delta \varphi _i^2>^n`$ $``$ $`\gamma ^{n1}<\delta \varphi _i^2>,`$ (A3)
$`<(\delta \varphi _i\delta \varphi _j)^2>^n`$ $``$ $`2^{n1}(\gamma \rho )^{n1}(<\delta \varphi _i^2>+<\delta \varphi _j^2>2<\delta \varphi _i\delta \varphi _j>).(ij)`$ (A4)
we get

$`{\displaystyle \frac{d\mu }{dt}}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n)}}{2^nn!}}\gamma ^n+w{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n)}}{n!}}(\gamma \rho )^n,`$ (A5)
which yields Eq. (11).

From Eqs. (8) and (9), we get

$`{\displaystyle \frac{d\delta \varphi _i}{dt}}`$ $`=`$ $`{\displaystyle \frac{d\varphi _i}{dt}}{\displaystyle \frac{d\mu }{dt}},`$ (A6)
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}F^{(2n+1)}{\displaystyle \frac{(\delta \varphi _i)^{2n+1}}{(2n+1)!}}+{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}F^{(2n)}[{\displaystyle \frac{\delta \varphi ^{2n}}{(2n)!}}{\displaystyle \frac{\gamma ^n}{2^nn!}}]`$ (A9)
$`+{\displaystyle \frac{w}{N}}{\displaystyle \underset{j}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}G^{(2n+1)}{\displaystyle \frac{(\delta \varphi _j\delta \varphi _i)^{2n+1}}{(2n+1)!}}`$
$`+{\displaystyle \frac{w}{N}}{\displaystyle \underset{j}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}G^{(2n)}[{\displaystyle \frac{(\delta \varphi _j\delta \varphi _i)^{2n}}{(2n)!}}{\displaystyle \frac{(\gamma \rho )^n}{n!}}]+\xi _i(t).`$
With the use of Eq. (A6), the calculation of $`d\gamma /dt`$ is performed as follows.

$`{\displaystyle \frac{d\gamma }{dt}}`$ $`=`$ $`{\displaystyle \frac{2}{N}}{\displaystyle \underset{i}{}}<\delta \varphi _i{\displaystyle \frac{d\delta \varphi _i}{dt}}>,`$ (A10)
$`=`$ $`{\displaystyle \frac{2}{N}}{\displaystyle \underset{i}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n+1)}}{(2n+1)!}}<\delta \varphi _i^{2n+2}>+{\displaystyle \frac{2}{N}}{\displaystyle \underset{i}{}}<\delta \varphi _i\xi _i>`$ (A12)
$`{\displaystyle \frac{2w}{N^2}}{\displaystyle \underset{i}{}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n+1)}}{(2n+1)!}}<\delta \varphi _i(\delta \varphi _i\delta \varphi _k)^{2n+1}>,`$
$`=`$ $`{\displaystyle \frac{2}{N}}{\displaystyle \underset{i}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n+1)}}{(2n+1)!}}B_{2n+2}<\delta \varphi _i^2>^{n+1}+2D`$ (A14)
$`{\displaystyle \frac{2w}{N^2}}{\displaystyle \underset{i}{}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n+1)}}{(2n+1)!}}B_{2n+2}<\delta \varphi _i(\delta \varphi _i\delta \varphi _k)>^{n+1}.`$
By using the mean-field approximation given by Eqs. (A2) and (A3) and

$$<\delta \varphi _i(\delta \varphi _i\delta \varphi _j)>^{n+1}(\gamma \rho )^n(<\delta \varphi _i^2><\delta \varphi _j>),(ij)$$
(A15)
we get

$`{\displaystyle \frac{d\gamma }{dt}}`$ $`=`$ $`2{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(2n+1)}}{2^nn!}}\gamma ^{n+1}2{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n+1)}}{n!}}(\gamma \rho )^{n+1}+2D,`$ (A16)
which leads to Eq. (12).

The calculation of $`d\rho /\underset{}{\text{t}}`$ is similarly performed by

$$\frac{d\rho }{dt}=\frac{1}{N^2}\underset{i}{}\underset{j}{}<\delta \varphi _i\frac{d\delta \varphi _j}{dt}+\frac{d\delta \varphi _i}{dt}\delta \varphi _j>,$$
(A17)
which yields Eq. (13).

For the AR model given by Eq. (14) with $`F(x)=1asin(x)`$ and $`G(x)=sin(x)`$, we get

$`F^{(\mathrm{})}(\mu )`$ $`=`$ $`casin(\mu ),\text{(}\mathrm{}=0\text{)}`$ (A18)
$`=`$ $`(1)^{n+1}asin(\mu ),\text{(}\mathrm{}=2n>0\text{)}`$ (A19)
$`=`$ $`(1)^{n+1}acos(\mu ),\text{(}\mathrm{}=2n+1\text{)}`$ (A20)
$`G^{(\mathrm{})}(0)`$ $`=`$ $`0,\text{(}\mathrm{}=2n\text{)}`$ (A21)
$`=`$ $`(1)^n,\text{(}\mathrm{}=2n+1\text{)}`$ (A22)
which yield Eqs. (15)-(17).

## B Derivation of Eqs. (54) and (55)

The moment method takes into account means, variances and covariances defined by

$`m_i`$ $`=`$ $`<\varphi _i>,`$ (B1)
$`C_{ij}`$ $`=`$ $`<\mathrm{\Delta }\varphi _i\mathrm{\Delta }\varphi _j>,\text{(}i,j=1\text{ to }N\text{)}`$ (B2)
where $`\mathrm{\Delta }\varphi _i=\varphi _im_i`$ and $`C_{ii}`$ denotes variances. By adopting the Gaussian decoupling approximation given by Eq. (10), we get

$`{\displaystyle \frac{dm_i}{dt}}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F_i^{(2n)}}{(2n)!}}B_{2n}<\mathrm{\Delta }\varphi _i^2>+{\displaystyle \frac{w}{N}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n)}}{(2n)!}}B_{2n}[<(\mathrm{\Delta }\varphi _k\mathrm{\Delta }\varphi _i)^2>]^n,`$ (B3)
$`{\displaystyle \frac{dC_{ij}}{dt}}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F_i^{(2n+1)}}{(2n+1)!}}B_{2n+2}(<\mathrm{\Delta }\varphi _j^2>^n+<\mathrm{\Delta }\varphi _i^2>^n)<\mathrm{\Delta }\varphi _i\mathrm{\Delta }\varphi _j>`$ (B6)
$`+{\displaystyle \frac{w}{N}}{\displaystyle \underset{k}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G^{(2n+1)}}{(2n+1)!}}B_{2n+2}[<(\mathrm{\Delta }\varphi _k\mathrm{\Delta }\varphi _j)^2>^n<\mathrm{\Delta }\varphi _i(\mathrm{\Delta }\varphi _k\mathrm{\Delta }\varphi _j)>`$
$`+<(\mathrm{\Delta }\varphi _k\mathrm{\Delta }\varphi _i)^2>^n<\mathrm{\Delta }\varphi _j(\mathrm{\Delta }\varphi _k\mathrm{\Delta }\varphi _i)>]+2D\delta _{ij},`$
where $`F_i^{(\mathrm{})}=F^{(\mathrm{})}(m_i)`$ and $`G_i^{(\mathrm{})}=G^{(\mathrm{})}(0)`$. By a proper re-arrangement, Eqs. (B3) and (B4) reduce to Eqs. (54)-(55).

| type of states | $`\zeta `$ | $`\delta \zeta `$ | $`\nu `$ | $`\sigma `$ |
| --- | --- | --- | --- | --- |
| S | F | 0 | 0 | F |
| R | 0 | 0 | F | 0 |
| P | F | F | F | F |

Table 1 Quantities in the S, R and P states of coupled AR networks: F and 0 denote finite and vanishing values, respectively.