Journal of the Physical Society of Japan

FULL PAPERS

Invariant Multiparameter Sensitivity to Characterize Dynamical
Systems on Complex Networks
Kenzaburo Fujiwara ∗ , Takuma Tanaka, and Kiyohiko Nakamura
Interdisciplinary Graduate School of Science and Engineering,Tokyo Institute of Technology, Yokohama 226-8502, JapanThe behavior of systems is determined by the parameters. Because we seldom know thedetailed structure of a system, metrics of parameter sensitivity should be independent of how
we model the system. We formulate a new parameter sensitivity metric, which we refer to as“invariant multiparameter sensitivity” (IMPS) because it gives the same result for a class ofequivalent models of a system. To investigate the property of IMPS, we firstly apply IMPS
to resistor circuits and linear dynamical systems. To examine the dependence of IMPS on
network structures, we secondly apply IMPS to nonlinear systems on complex networks. We
find that the IMPS of networks of phase oscillators is essentially independent of the numberof oscillators. We examine the network-structure dependence of IMPS using a simplified
solvable model.

1. Introduction
Over the past few decades, research has increasingly focused on social, biological, economical, ecological, and geological systems.1–5) The outputs of these systems depend on parameters. In these systems, the parameters may fluctuate stochastically and undergo long-term
drift. Thus, understanding the effect of variations in the parameters on the system outputs isvital for understanding their dynamics.6, 7) For example, this approach is used to characterize
electric circuits that are insensitive to variations in individual circuit elements.7–9)
In this study, we focus on the changes in the response of system outputs when the parameters are changed; this relation is referred to as the “parameter sensitivity” of the system.
Parameter sensitivity has been intensively studied in circuit theory, particularly in resistor–
capacitor networks.7–10) Although parameter sensitivity was originally proposed to evaluate
the sensitivity of electric circuits, it has also been used in many other fields, such as biology,11–13) ecology,14) and economics.15) Previous studies have proposed several metrics for∗kenzaburo.fujiwara@gmail.com
1/13

J. Phys. Soc. Jpn.

A

H

FULL PAPERS

B

C
R

unknown system

E

Fig. 1.

R’1

R’m(m− 2) resistors

E

E(Color online) Unknown heating system (A) and series resistor models (B,C). The models B and C

are electrically equivalent when R′1 + R′2 + · · · + R′m = R.

parameter sensitivity.16, 17) Single parameter sensitivity (SPS) is the ratio of the change in output to the change in parameter.18) Since most systems have more than one parameter, generalizing SPS to multiple parameters is the next logical step. Multiparameter sensitivity (MPS)
is a generalization of SPS for multiple parameters.18) MPS is defined as the square root of thesum of the squares of SPSs of the system. MPS has been used in analyzing networks withcircuit elements.9, 18, 19) In biology, MPS was used in clarifying how the integration of molecular components generates robust systems such as the Escherichia coli heat-shock response
system and the Drosophila melanogaster circadian clock.11)
The parameter sensitivity of a system is often numerically estimated by quantifying theparameter sensitivity of its models.8, 11–13) Because, in most cases, we do not know the structures and dynamics of systems completely, a wide variety of models are proposed for a single
system. As an example, let us assume that we want to quantify the parameter sensitivity of thesystem of Fig. 1A. When connected to a voltage source E, this system exhibits heat flow H.
Here, we assume that this system comprises linear resistors. Even under this assumption, thissystem can be modeled by a large number of electric circuit models, such as those in Figs. 1Band 1 C. Unfortunately, as will be shown in the next section, MPS gives different results forthe models in Figs. 1B and 1 C, which are electrically equivalent models of a single system ofFig. 1A. Since MPS depends on the choice of the models, we cannot quantify the parameter
sensitivity of the system using MPS.
To address this problem, we propose a new sensitivity metric, which we call “invariant
multiparameter sensitivity” (IMPS). In Sect. 2, we briefly review previously proposed sensitivity metrics and then define IMPS. We use the example of serial resistor models to present“invariance,” which is an important property of IMPS. In Sect. 3, IMPS for linear dynamical systems is examined. Furthermore, we study the IMPS for phase oscillator models oncomplex networks in Sect. 4. IMPS is suggested to vary according to network structures. Finally, in Sect. 5, we summarize the main results and discuss the potential applications andlimitations of IMPS. Detailed derivations can be found in Supplemental Materials.20)
2/13

J. Phys. Soc. Jpn.

FULL PAPERS

2. Invariant Multiparameter Sensitivity
Models of dynamical systems are expressed by the first-order differential equations~x˙ = F t, ~x, ~p ,(1)

where t is the time, ~x = [x1 , x2 , x3 , . . . , xn ] is the state variable vector, and ~p =[p1 , p2 , p3 , . . . , pm ] is the parameter vector. Let q be the output whose sensitivity to parameters we want to quantify. Here, the output can be any value derived from the system, such as
oscillation period, cell growth rate, or energy consumption. For a given output q, the single
parameter sensitivity with respect to the parameter pi , which we call SPSi , is defined as∂ ln q
pi ∂q
=
.
(2)
SPSi =q ∂pi ∂ ln pi
SPSi is the ratio of the change in the output q to the change in the parameter pi .18) SPSi does
not quantify the sensitivity of q with respect to any other parameters. Generalizing SPS tomultiple parameters, Goldstein and Kuo defined MPS as the L2 norm of SPSs:18)
m
X
2
MPS =
SPS2i .

(3)i=1

MPS has been used to quantify the sensitivity of the output to changes in the entire parameterset of models.9, 18, 19) However, MPS often gives different results for two equivalent models ofa system. Consider a model of a circuit with one resistor R, as shown in Fig. 1B. Let H denotethe heat flow from the circuit and let E denote the voltage of the voltage source. We assume
that the heat flow H = E 2 /R is the output of this model. Because the SPS of this circuit withrespect to R is unity, the MPS is also unity. A current-voltage relationship identical to thatof Fig. 1B can also be implemented by equivalent circuits, such as that shown in Fig. 1C(assuming R′1 + R′2 + . . . + R′m = R). The MPS of the circuit of Fig. 1C is given by!
mX
R′i ∂H 22
MPS =H ∂R′i
i=1
=

m
X
R′ 2
i

R2
i=1

.

(4)

Hence, the MPS of this model decreases as m increases.
In contrast, the L1 norm of SPSs of the circuit for Fig. 1C is given bymm
X
X
R′i ∂H|SPSi | =H ∂R′i
i=1
i=1
= 1,

(5)

which equals the L1 norm of SPSs for the circuit of Fig. 1B. Thus, we introduce a new metric

3/13

J. Phys. Soc. Jpn.

FULL PAPERS

IMPS, which is the sum of all absolute values of SPSs:m
X
IMPS =|SPSi | .

(6)

i=1

The IMPS of the circuit is constant, regardless of m because the output H is a homogeneous function of degree −1 of resistances R1 , R2 , . . ., and Rm . Indeed, assuming that theoutput q(p1 , p2 , p3 , . . . , pm ) is a homogeneous function of the degree k and that the SPSs in
Eq. (6) have the same sign, we obtain
IMPS =

m
X
pi ∂q
i=1

=

q ∂pi

|k| ,

where we use Euler’s theorem!
mX∂q= kq(p1 , p2 , p3 , . . . , pm ).
pi
∂p
i
i=1

(7)

(8)

This is also the case with equivalent circuits of resistors in parallel. These results suggestthat IMPS is independent of the models, and may be a metric for quantifying the sensitivityof the system itself rather than the model. The IMPSs are the same for all models satisfying
the following conditions: (i) the output is a homogeneous function of the parameters and (ii)the SPSs have the same sign. These conditions are not necessarily satisfied. However, in thefollowing, we show that these conditions are approximately satisfied in several models of
linear and nonlinear dynamics.
3. Linear Dynamics
There are many systems that can be modeled by elements interacting with each other.21–24)
As a simple example of these models, we examine a model consisting of N particles. Assume
that particle i has the position xi and the velocity ẋi , and is attracted to other particles via theconnection weight matrix W = (Wi j ). The matrix W is symmetric, has no negative elements,
and is assumed to define a connected graph of particles. Then, the dynamics are defined byNXdxi= bi +
Wi j (x j − xi ),
ẋi =
dt
j=1

(9)

where bi is the velocity of the particle i with no interaction. The velocity bi is drawn from aunit normal distribution. The positions of the particles converge to fixed points. The variance
V of the positions at the steady state is used as the output of the model, whose sensitivity is
to be measured. We assume a connected graph because the variance cannot be defined whenparticles become disconnected.

4/13

J. Phys. Soc. Jpn.

FULL PAPERS

B 2.1

0.2

IMPS
MPS
1

0

20

40

60

80

100

Number of Particles N

Fig. 2.

2

1.9
1000

0.1

3000

MPS

IMPS
MPS

2

IMPS

IMPS & MPS

A

0
5000

Number of Connections

IMPS (solid line) and MPS (dashed line) for linear dynamical models. Panel A shows IMPS and MPS

as functions of N. Panel B shows IMPS and MPS as functions of the number of connections for N = 100.

By using the Laplacian matrix L of a weighted graph25) defined by the connection weightmatrix W, we obtain
IMPS =!
1 X Wi j ~ T −1 ∂L̃ −1−1 ∂L̃
−L̃ + L̃L̃−1~bb L̃
N hi ji
V
∂Wi j
∂Wi j

1 β ~ T −1
− b L̃ (1L̃−1 + L̃−1 1)L̃−1~b ,(10)
N V
where hi indicates the summation over the connected particle pairs, and L̃i j = Li j + β. β , 0 is+

introduced so that L̃ is invertible. IMPS is independent of β. The detailed derivation is given
in Supplemental Materials Analysis S1.20) By assuming that all SPSs have the same sign, we
obtain IMPS = 2 because the variance is a homogeneous function of degree −2, regardless of
network size and the structure given by the connection weight matrix W as shown in Eq. (S20).
Let us examine whether the result IMPS = 2 holds in numerical simulations. Figure 2A
shows that this relation holds for the network of all-to-all coupling, Wi j = 0 for i = j and
Wi j = 1 otherwise, for a wide range of N. This suggests that all SPSs have the same sign inthis model. In contrast, the MPS decreases as N increases (Fig. 2A).
In the simulation of Fig. 2B, we start from a fully connected system of size N = 100
and eliminate a connection at each step. We confirmed that W constituted a connected graph.
Thus, the graph at each step is an Erdős–Rényi random graph.26) The IMPS and MPS are
shown in Fig. 2B as functions of the number of connections. This figure shows that MPS depends much more strongly on the number of connections than IMPS. Therefore, these results
indicate that the condition that all SPSs have the same sign is almost satisfied. Moreover,these results indicate that, for these linear models, the IMPS quantifies how sensitively thesystem’s output, that is, the variance of the particle positions, changes in response to a parameter change in a manner that is independent of N and the number of connections.
5/13

J. Phys. Soc. Jpn.

FULL PAPERS

4. Phase Oscillator Networks
Let us consider a system of oscillatory elements on a network. Oscillatory elements areoften modeled by phase oscillators.21, 27, 28) Because the precise structure of a network may be
unknown to us, to investigate how IMPS is affected by the choice of network structure used tomodel the system, we compare IMPSs of phase oscillators on the three types of network, thatis, Barabási–Albert, regular random, and Watts–Strogatz networks. Now, consider a model of
N oscillators connected to each other by the adjacency matrix A = (Ai j ), where Ai j = 1 forconnected pairs and Ai j = 0 otherwise. The dynamics of the oscillator i are described byNXdθi= ωi +
Ki j sin(θ j − θi ),
dt
j=1

(11)

where ωi is the natural frequency of the oscillator i and Ki j is the element (i, j) of the connection weight matrix K.21, 29) For the numerical simulation, we use K = αA, where α > 0 is thecoupling strength, and draw the natural frequency ωi from a unit normal distribution. Without
PN
loss of generality, we can assume that i=1
ωi = 0. For α, we use a sufficiently large value sothat all phase oscillators are phase locked. We use the circular variance Vc of the oscillatorsin the phase-locked state30) as the output:
1√ 2
C + S 2,
Vc = 1 − r = 1 −N
PNPNsin θi .
cos θi , and S = i=1
where r is the Kuramoto order parameter, C = i=1(12)

To simulate phase oscillators on complex networks, we use the adjacency matrices given

by Barabási–Albert31) (a well-known model of scale-free networks) and regular random networks with an average degree of 4. We use output Vc at t = 100 for Barabási–Albert andregular random networks. To generate a Barabási–Albert network, we use Gephi complex
generators. We start from 2 connected vertices and add a vertex with 2 edges in each step until we have N vertices. We generate regular random networks with degree 4 using the Matlab
Random Regular generator. The phases of all oscillators are distributed uniformly at t = 0.
By using the implicit function theorem, we derive ∂Vc /∂Ki j [Eq. (S-32)].20) We calculate
the IMPS numerically using Eq. (S-33).20) Unlike for equivalent electric circuits and linear
dynamical models, IMPSs for phase oscillator networks are not constant (Fig. 3). This result
is attributed to Vc not being a homogeneous function. Figure 3A shows the results of 10 000
phase oscillators on 5 Barabási–Albert networks and Fig. 3B shows the results of 10 000
phase oscillators on 5 regular random networks generated by different random seeds. The
IMPSs converge to 2, as shown in Figs. 3A and 3B. This is because, as the coupling strength
α increases, the models can be described by linearized dynamics. The nonlinear dynamics of

6/13

J. Phys. Soc. Jpn.

FULL PAPERS
Table I.

Mean
Standard Deviation

Average path length of oscillator networks (N = 10 000)

Barabási–

Regular

300-connection

3000-connection

AlbertRandomRewired Watts–Strogatz

Rewired Watts–Strogatz

36.1

10.9

0.343

2.88 × 10−2

5.03

7.73
−2

2.20 × 10−3

1.53 × 10

Eq. (11) become the linear dynamics of Eq. (9), and the circular variance of Eq. (12) becomes
proportional to the variance in the large α limit.
In Fig. 4, for models on Barabási–Albert networks, we compare the IMPS and MPS as
functions of coupling strength α for models with different N values. As shown in Fig. 4A,
the IMPSs for models with different N values are close together. In contrast, as shown in
Fig. 4B, the MPSs for models with different N values are well separated. IMPSs for models
with different N values converge to the same value, 2, for the nonlinear models as well as thelinear models.
Next, to examine another type of complex network, we use the adjacency matrices given
by Watts–Strogatz networks1) with an average degree of 4. We use the output Vc at t = 200
for Watts–Strogatz networks. To generate a Watts–Strogatz network,1) we start from a regularring lattice with N vertices; each vertex is connected to 4 neighbors, 2 on each side. We
randomly select an edge and rewire one of its endpoints step by step either 300 or 3000
times. The phases of all oscillators at t = 0 are drawn from a Gaussian distribution with amean of 0 and a standard deviation of 0.1.
Figure 3C1 (3C2) shows the IMPS for 10 000 phase oscillators on Watts–Strogatz networks with 300 (3000) of 20 000 connections rewired. Figures 3C1 and 3C2 are generatedusing different random seeds. For the models on Watts–Strogatz networks, although the IMPS
approaches 2 as α increases, the IMPS deviates from 2 more than the IMPS for the models on
Barabási–Albert and regular random networks. As shown in Fig. 3C1, the IMPS curves are
distributed more widely than those of Figs. 3A and 3B. The result suggests that the IMPS ofthe models on Watts–Strogatz networks varies from network to network. On the other hand,
IMPSs are within a narrower range for Barabási–Albert and regular random networks. For
the models on Barabási–Albert and regular random networks, the IMPS converges to 2 fromabove, whereas, for models on Watts–Strogatz networks, the IMPS converges to 2 from below. These observations suggest that the circular variance for Barabási–Albert and regularrandom networks changes much more sensitively depending on the coupling strength of each
7/13

J. Phys. Soc. Jpn.

FULL PAPERS

A

B
2.2

IMPS

IMPS

2.2
2
1.8

1.8

1.6

1.6

2

3

4

5

α

6

7

2

8

3

4

5

6

7

8

α

C2

C1
2.2

2.2

IMPS

IMPS

2

2
1.8

2
1.8
1.6

1.6

5

40 45 50 55 60 65 70

10 15 20 25 30 35

αFig. 3.

α

IMPS for 10 000 oscillators on complex networks. IMPSs for the models on Barabási–Albert (A),regular random (B), and 300-connection rewired Watts–Strogatz (C1), and 3000-connection rewired Watts–
Strogatz (C2) networks are shown. In panels A and B, the coupling strength α varies from 2 to 8 in steps of 0.2.
In panel C1, α varies from 40 to 70 in steps of 1. In panel C2, α varies from 5 to 35 in steps of 1.

A

B
0.4
N = 100
N = 1000
N = 10000

N = 100
N = 1000
N = 10000

0.3

MPS

IMPS

2.2

2.1

0.2
0.1

2

2

4

6

α

Fig. 4.

8

10

0

2

4

6

α

8

10

IMPS and MPS for 100, 1000, and 10 000 oscillators on Barabási–Albert networks. (A) IMPS and(B) MPS as functions of coupling strength α for the models on Barabási–Albert networks. In both panels, thecoupling strength α varies from 2 to 10 in steps of 0.2.

vertex than the circular variance for Watts–Strogatz networks when the coupling strength α
is small. This implies that IMPS reflects some difference in network structures and network
metrics among Barabási–Albert and regular random networks and Watts–Strogatz networks.
To clarify the factor reflected in IMPS, next we examine a solvable oscillator model. Compared with Barabási–Albert and regular random networks, Watts–Strogatz networks have
longer average path lengths (Table I). Watts–Strogatz networks have a larger number of long
one-dimensional (1-D) lattice like structures (Fig. 5A) than Barabási–Albert and regular random networks. To investigate the effects of 1-D lattice like structures, we simplify a regular
1-D lattice to a path graph, for which all vertices and edges lie on a single straight line
8/13

J. Phys. Soc. Jpn.

FULL PAPERSB
ceructure
e St
k
i
L

1
2

S

1-D
La
tti

A

ut
rtc
ho
N −1
N

C
2.2

N=2
N=3
N=4
N=5

2

IMPS

N=6
N=7

1.8
N=8
N=9

1.6

N = 10

1.4
0

5

10

α

Fig. 5.

(Color online) One-dimensional lattice like structure, path graph, and IMPS for path graphs of different

lengths. Panel A shows a schematic of a 1-D lattice like structure. Panel B shows a path graph with N vertices.
Panel C shows the IMPS for models with different-length path graphs as a function of the coupling strength α.
In panel C, the phases of N oscillators are separated by regular intervals.

(Fig. 5B). When N phase oscillators are connected to each other along a path graph and havenatural frequencies −ω1 = ωN = 1 and ωi = 0 for 1 < i < N, their phases are separated byregular intervals (∆θ). IMPSs for the oscillators on different lengths can be obtained as! zN − 1 −zN + N(z − 1) + 1 tan(∆θ) ,(13)
Im IMPS = 2
N r(1 − r)z−1
(z − 1)2where r is the Kuramoto order parameter and z = ei∆θ (see Supplemental Materials Analysis
S3 for derivation20) ). As shown in Fig. 5C, the length of the path graph determines whether the
IMPS converges from above or below. The models of shorter (longer) path graphs correspondto models on Barabási–Albert and regular random networks (Watts–Strogatz networks). To
further examine how the dependence of IMPS on α is affected by the average path length,
we plot the difference between IMPS for a small α and that for a large α (Fig. S1).20) As
is consistent with Fig. 5C, the correlation coefficient r between the average path length and
9/13

J. Phys. Soc. Jpn.

FULL PAPERS

IMPS(α = 10) − IMPS(α = 100) on 50 Watts–Strogatz networks (N = 1000, rewiring

probability= 0.05) is negative (r = −0.43, p < 0.01). Other network metrics such as diameterand radius are also correlated to IMPS’s difference. On the other hand, for Barabási–Albert
networks, the average path length is not significantly correlated to IMPS(α = 2.5)−IMPS(α =100), probably owing to the narrow range of the average path length of randomly generated
50 Barabási–Albert networks (N = 1000). These observations suggest that the networks
with shorter average path lengths make the model more sensitive to changes in the couplingstrength of each vertex.
5. Discussion
In this paper, we reviewed the previously proposed metrics SPS and MPS, which quantify the sensitivity of the output to the values of parameters. We formulated an improved
metric called IMPS, which we applied to series-resistor circuits and linear dynamical systems. By calculating the IMPS for series-resistor circuits, we found that IMPS gives identical results for equivalent circuits of Figs. 1B and 1C. The analysis and simulation of linear
dynamical models revealed that IMPS is less dependent on the number of particles in thelinear dynamical model. We also calculated IMPS for nonlinearly coupled phase oscillators on Barabási–Albert, regular random, and Watts–Strogatz networks. In all of these cases,
the IMPS converges to values close to 2. In the limit of a large coupling strength, IMPS isless dependent on the type of network, that is, the Barabási–Albert network, regular random
network, Watts–Strogatz network, and path graph. However, for a small coupling strength,
IMPSs differ among the types of network used to model the unknown system. Moreover, fora small coupling strength, average path length negatively correlates to IMPS. Thus, for systems of phase oscillators on networks, our results show that the system with a shorter averagepath length is more sensitive to parameter changes.
The IMPS gives identical results for equivalent circuits (Fig. 1). The invariance of IMPS
was previously reported only for resistor–capacitor networks.7–10) The analysis of the circuits
indicates that dividing an element into several elements or combining several elements into asingle element does not change the IMPS. This is also the case with the linear and nonlinear
dynamical systems. These results suggest that IMPS can be a more appropriate metric than
MPS in quantifying the sensitivity of unknown systems, because the precise number of elements in the system in question is often unknown. However, the conditions for the invarianceof IMPS may be satisfied only in a limited class of systems.
To use IMPS to compare the sensitivity of various systems, two points must be considered:

10/13

J. Phys. Soc. Jpn.

FULL PAPERS(i) IMPS depends on the shift of the zeroes of parameters and (ii) IMPS depends on thedimension of the units of parameters. As an example of the first point, assume the pressure P
as the output of a system consisting of an ideal gas in a cuboid. The ideal gas formula is
nRT
P=
,
(14)
V
where V is the volume of the gas, n is the number of moles of the gas in the cuboid, R is
the ideal gas constant, and T [K] is the absolute temperature in Kelvin of the gas. Here, we
assume that R is a constant and n, T , and V are the parameters. The IMPS for this model is
3. If we give the temperature C in Celsius, we have T = C + 273. The IMPS for this modelC|. As an example of the second point, denote the width, depth,
expressed in Celsius is 2+| C+273

and height of the cuboid by w, d, and h, respectively. This gives
nRT
.
(15)
P=
wdh
The IMPS for this model expressed in terms of n, T , w, d, and h is 5. Thus, the points ofzeroes of the parameters and the dimensions of units should be consistent between the models
describing systems. SPS and MPS also suffer from the same difficulty.
To summarize, in this paper, we investigated the IMPS with respect to the output of linear
and nonlinear models that can be approximately linearized. In strongly nonlinear systems,
such as chaotic systems, how their nonlinearity is reflected in IMPS should be examined. In
addition, the relationship between IMPS and the Lyapunov exponent is of interest. Future
work should also consider the IMPS of other types of nonlinear system.
Acknowledgment
This work was supported by MEXT/JSPS KAKENHI Grant Numbers 24651184 and25115710. A preliminary version of this paper was presented at ICONIP 2014.

11/13

J. Phys. Soc. Jpn.

FULL PAPERS

References
1) D. J. Watts and S. H. Strogatz: Nature 393 (1998) 440.
2) R. V. Solé and J. M. Montoya: Proc. R. Soc. Lond. B 268 (2001) 2039.
3) B. S. Frey and A. Stutzer: J. Econ. Lit. 40 (2002) 402.
4) H. Ukai, T. J. Kobayashi, M. Nagano, K. Masumoto, M. Sujino, T. Kondo, K. Yagita,
Y. Shigeyoshi, and H. R. Ueda: Nat. Cell biol. 9 (2007) 1327.
5) H. Cheng, R. L. Edwards, W. S. Broecker, G. H. Denton, X. Kong, Y. Wang, R. Zhang,and X. Wang: Science 326 (2009) 248.
6) M. E. Csete and J. C. Doyle: Science 295 (2002) 1664.
7) J. V. Leeds and G. Ugron: IEEE Trans. Circuit Theory 14 (1967) 188.
8) P. Goddard and R. Spence: Electron. Lett. 5 (1969) 351.
9) A. Rosenblum and M. Ghausi: IEEE Trans. Circuit Theory 18 (1971) 592.
10) T. Roska: Electron. Lett. 4 (1968) 281.
11) K. Maeda and H. Kurata: J. Ther. Biol. 272 (2011) 174.
12) H. Kurata, T. Tanaka, and F. Ohnishi: PLoS ONE 2 (2007) e1103.
13) J. Yoon and T. S. Deisboeck: PLoS ONE 4 (2009) e4560.
14) K.-E. Lindenschmidt: Ecol. Model. 190 (2006) 72.
15) P. Vermeulen and D. De Jongh: Appl. Math. Model. 1 (1976) 29.
16) J. C. Helton: Reliab. Eng. Syst. Safe. 42 (1993) 327.
17) I. M. Sobol’: Math. Comput. Simul. 55 (2001) 271.
18) A. Goldstein and F. Kuo: IRE Trans. Circuit Theory 8 (1961) 177.
19) J. D. Schoeffler: IEEE Trans. Circuit Theory 11 (1964) 271.
20) (Supplemental Materials) Analyses for linear dynamical models and nonlinear models
are provided online. Analyses for phase oscillators on path graphs and a figure for thecorrelation between IMPS and the average path length are also provided online.
21) Y. Kuramoto: Chemical Oscillations, Waves and Turbulence (Springer, Berlin, 1984).
22) S. Nakata, Y. Doi, and H. Kitahata: J. Phys. Chem. B 109 (2005) 1798.
23) H. Kitahata, J. Taguchi, M. Nagayama, T. Sakurai, Y. Ikura, A. Osa, Y. Sumino,
M. Tanaka, E. Yokoyama, and H. Miike: J. Phys. Chem. A 113 (2009) 8164.
12/13

J. Phys. Soc. Jpn.

FULL PAPERS

24) M. Yomosa, T. Mizuguchi, and Y. Hayakawa: PloS ONE 8 (2013) e81754.
25) R. Olfati-Saber and R. M. Murray: IEEE Trans. Automatic Control 49 (2004) 1520.
26) P. Erdős and A. Rényi: Publ. Math. Debrecen 6 (1959) 290.
27) A. Takamatsu, T. Fujii, and I. Endo: Phys. rev. lett. 85 (2000) 2026.
28) I. Aihara, H. Kitahata, K. Yoshikawa, and K. Aihara: Artif. Life Robotics 12 (2008) 29.
29) T. Ichinomiya: Phys. Rev. E 70 (2004) 026116.
30) N. I. Fisher: Statistical Analysis of Circular Data (Cambridge University Press, 1995).
31) A.-L. Barabási and R. Albert: Science 286 (1999) 509.

13/13