Detecting connectivity of small, dense oscillator networks from dynamical measurements based on a phase modeling approach

JID:PLA AID:21972 /SCO Doctopic: Nonlinear science

[m5Gv1.5; v 1.96; Prn:15/05/2013; 16:08] P.1 (1-6)

Physics Letters A ••• (••••) •••–•••

1

Contents lists available at SciVerse ScienceDirect

2

67 68

3

Physics Letters A

4 5

69 70 71

6

72

www.elsevier.com/locate/pla

7

73

8

74

9

75

10

76

11 12 13

Detecting connectivity of small, dense oscillator networks from dynamical measurements based on a phase modeling approach

14 15 16 17 18

Isao T. Tokuda

b

, Mahesh Wickramasinghe , István Z. Kiss

a

a

Department of Mechanical Engineering, Ritsumeikan University, 1-1-1 Noji-higashi, Kusatsu, Shiga 525-8577, Japan b Department of Chemistry, Saint Louis University, 3501 Laclede Ave., St. Louis, MO 63103, USA

23 24 25 26 27 28 29 30 31 32

79 81 82 83 84 85

20 22

78 80

a,∗

19 21

77

86

a r t i c l e

i n f o

Article history: Received 14 February 2013 Received in revised form 6 May 2013 Accepted 8 May 2013 Available online xxxx Communicated by C.R. Doering Keywords: Oscillator networks Connectivity Phase equation Data analysis

a b s t r a c t

87 88

An approach is presented for detecting the connectivity between the oscillator elements from the measured multivariate time series data. Our methodology is based upon the phase equation modeling of the oscillator networks, where not only the connection matrix but also the natural frequencies and the interaction function of the oscillators are estimated. Application of this technique to simulated data as well as experimental ones from electrochemical oscillators shows its capability for precise detection of defects in the connection matrix for small-size networks. Dependence of the methodology on the observational noise, the network size, the number of defects, and the data length is also examined. © 2013 Published by Elsevier B.V.

89 90 91 92 93 94 95 96 97 98

33

99

34

100

35 36

1. Introduction

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

A network of interacting oscillators can be found in diverse fields of natural science and engineering [1–3]. Here, connectivity of the network elements plays a crucial role on the formation of the collective dynamics. Effect of the complex network topology on the dynamics of coupled oscillators has been intensively investigated [4]. To deal with the real-world systems, however, it is extremely rare that the detailed connectivity of the network can be directly and noninvasively investigated. In many systems especially in biology, direct measurement of the connectivity may destroy the true nature of the coupling functions that underlie the real network dynamics. For instance, in the studies of circadian rhythms, fundamental mechanisms for connecting the circadian cells in the superchiasmitic nucleus remain largely unknown [5,6]. Detailed physiological investigation of the in vitro tissue may not reveal the true network function that is inherent in situ. Because of these difficulties, estimation of the network connectivity from time series data recorded in a noninvasive manner is an awaited technique. Along this line, several approaches have been proposed up to date. Information transfer [7], mutual predictability [8], recurrence properties [9], and permutation-based asymmetric association measure [10] have been utilized to identify the coupling directions. Index for partial phase synchronization has been developed to distinguish direct from indirect interactions [11–13].

61 62 63 64 65 66

*

Corresponding author. Tel.: +81 77 561 2832; fax: +81 77 561 2665. E-mail address: [email protected] (I.T. Tokuda).

0375-9601/$ – see front matter © 2013 Published by Elsevier B.V. http://dx.doi.org/10.1016/j.physleta.2013.05.016

A graph theoretic approach to detect causalities in multivariate time series has been recently developed [14]. Response properties of the network dynamics to external stimuli have been also exploited [15,16]. Here, we focus on a phase description of the oscillator networks [2]. Under certain conditions, which are described in detail in Section 2, a network of weakly coupled self-sustained oscillators can be reduced to a dynamics of interacting phase oscillators. This drastically simplifies the modeling assumptions and enables straightforward estimation of the phase dynamics. We have applied the multiple-shooting method to fit the phase equations to multivariate time series to show that the interaction functions as well as the natural frequencies can be well identified from a globally coupled populations [17]. Kralemann et al. [18,19], on the other hand, carried out a data fitting using the probability density function of a modified phase to estimate the phase dynamics from time series and Blaha et al. [20] applied it to two interacting electrochemical oscillators. Our approach has the practical advantage of algorithmic simplicity and hence it is straightforward to implement. It has been, however, applied only to the simplest case of all-to-all coupling, where all oscillator elements are connected to all the others. The aim of the present Letter is to extend our previous approach to a network dynamics of an arbitrary topology. In particular, we focus on the case that the network has zero-or-unity connections, where the zero connections are referred to as defects in the coupling. We examine the capability of our approach to reveal the network connectivity from multivariate data recorded from coupled systems.

101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:PLA

AID:21972 /SCO Doctopic: Nonlinear science

1

[m5Gv1.5; v 1.96; Prn:15/05/2013; 16:08] P.2 (1-6)

I.T. Tokuda et al. / Physics Letters A ••• (••••) •••–•••

2

2. Method

2 3 4

Our problem can be stated as follows. Consider a system of N weakly coupled nearly identical limit cycle oscillators:

5 6 7

x˙ i = F i (xi ) +

C N

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

N 

T i , j G (xi , x j ),

(1)

j =1, j =i

where xi and F i (i = 1, 2, . . . , N) stand for state variables and autonomous dynamics of the ith oscillator, respectively, C represents the coupling constant, and G is an interaction function between ith and jth oscillators. The matrix { T i , j } describes connectivity between the oscillators, where the present problem is applicable to both uni- and bi-directionally coupled oscillators. Our assumption is that in isolated condition (i.e., C = 0) each oscillator F i generates a stable limit cycle with similar natural frequencies ωi . The functions F i should be similar among the N oscillators near the limit cycle trajectory. Then the phase reduction theory [2] states that for a weak coupling C the network  N dynamics can be reduced to C the phase equations: θ˙i = ωi + N j =1, j =i T i , j H (θ j − θi ) (θi : phase of ith oscillator; H : interaction function). As a recording condition, we assume that simultaneous measurement of all oscillators is made as {ξi (nt ) = g (xi (nt )): n = 1, . . . , M }iN=1 (g: observation function, t: sampling time). Our objective of inferring the network connectivity through recovering the phase dynamics is accomplished under conditions that (i) the underlying dynamics (1) are unknown, (ii) the interaction function is based upon a difference coupling (i.e., H (0) = 0), (iii) the coupling constant C associated with the measured data is in a non-synchronized regime, and (iv) the connection matrix is composed of zero-or-unity elements (i.e., T i , j = 0 or 1), where most of the connections exist (non-sparse matrix). Number of the defects (T i , j = 0) is denoted by β . Concerning the condition (ii), the difference coupling provides a good assumption, because, in many systems, the phase interaction disappears when the oscillators are in a complete in-phase relationship. Concerning the condition (iv), such a non-sparsely connected system can be found in real-world systems. For instance, all-toall connections are assumed in various systems including arrays of Josephson junctions [21] and lasers [22], population of circadian oscillators [5,6], ensembles of electrochemical oscillators [23], and neuronal populations [24]. It can easily happen that a small portion of the connections is destroyed by some damage to the system, resulting in a non-sparsely connected system. Our approach to the problem is based upon the following four steps.

data. Simultaneous estimation of the unknown parameters all at once is, however, an ill-conditioned problem because of the redundancy of the unknown parameters. Therefore, we divide the ˜ i , a j , b j } and { T˜ i , j } and estiparameters into two groups as {ω mate them separately in the next two steps. 3. In the first step, the connection matrix is assumed to be of all-to-all type ( T˜ i , j = 1) and estimate the rest of the param˜ i , a j , b j }. This provides a reasonable eters denoted by p = {ω assumption for densely connected systems, in which most of the oscillators are coupled to each other. The parameters p are estimated by the multiple-shooting method [25]. We denote the time evolution of the phase equations (2) with respect to an initial condition θ by φ t (θ , p ). Then, at each sampling time t = nt, the phase equation must satisfy the boundary conditions: θ((n + 1)t ) = φ t (θ(nt ), p ). With respect to the unknown parameters p, we solve these nonlinear equations by the generalized Newton method. The evolution function φ t is integrated numerically. For the computation of the gradients ∂φ/∂ p which are needed for the Newton method, variational equations of the phase equations (2) are also solved numerically. A necessary condition to solve the nonlinear equations is that the number of the unknown parameters p is less than the number of the equations, that is, N + 2D < N ( M − 1). This always holds in the case we have enough data points M. ˜ i , a j , b j } estimated in 4. In the second step, the parameters {ω the previous step are fixed. With respect to the connection matrix p = { T˜ i , j } as the rest of the unknown parameters, the multiple-shooting method was applied again in a similar manner as in the previous step.

50 51 52 53

1. Extract phase signals θi (t ) from the data ξi (t ). The phase θ is defined by a simple piecewise linear formula, in which it is increased by 2π at every local maximum of ξ(t ) and between the local maxima it grows in proportion to time [3]. 2. Fit the phases {θi (t )} to the phase equations:

54 55 56 57 58 59 60 61 62 63 64 65 66

θ˙i = ω˜ i +

C N

N 

˜ (θ j − θi ), T˜ i , j H

(2)

j =1

˜ i } and { T˜ i , j } represent approximate values for the where {ω natural frequencies and those for the connection matrix, re˜ stands for an approximate function for the interspectively. H action H , which is in general nonlinear and periodic with respect to 2π . The approximation is basedon a Fourier ex˜ (θ) = Dj=1 [a j sin j θ + pansion up to the order of D as H b j (cos j θ − 1)]. The above phase equations have unknown parameters of {ω˜ i , a j , b j , T˜ i , j }, which should be optimized to be fitted to the

68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

It should be noted that, in the above method, the first and the second steps can be repeated in an iterative manner to improve the estimates. We do not employ such a procedure, because our preliminary study indicated that the iterative algorithm simply increased the computational cost but did not show a clear improvement.

98 99 100 101 102 103 104

3. Results

105 106

We applied this technique to a prototypical example of weakly coupled limit cycle oscillators. We considered the following system of Rössler equations with diffusive coupling:

y˙ i = αi xi + 0.15 y i + z˙ i = 0.2 + zi (xi − 2),

107 108 109 110

x˙ i = −αi y i − zi ,

48 49

67

111

C N

N 

112

T i , j ( y j − y i ),

113 114

j =1

(3)

where i = 1, . . . , N. To consider an inhomogeneity of the network elements, parameter values αi , which determine rotation speed in the (x, y )-space, were varied among the oscillators as αi = 1 + 0.01 · i (i = 1, . . . , N). Each Rössler oscillator generates a limit cycle attractor for the chosen parameter values in the absence of coupling (i.e., C = 0). The multivariate data were recorded as { y i (t )}iN=1 . We started with the case of N = 5. As the coupling matrix, all-to-all connections were assumed (T i , j = 1). The data { y i (t )}5i =1 were recorded at a coupling strength of C = 0.01, which corresponds to non-synchronized regime. The sampling interval was set to be t = 0.32 for the extraction of the phases {θi (t )}. Then to apply the multiple-shooting method the data have been down sampled to t = 250 · 0.32 and the total of 1000 data points have been collected for the parameter estimation. As an initial condition, the unknown parameter values were all set to be zero, i.e.,

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:PLA AID:21972 /SCO Doctopic: Nonlinear science

[m5Gv1.5; v 1.96; Prn:15/05/2013; 16:08] P.3 (1-6)

I.T. Tokuda et al. / Physics Letters A ••• (••••) •••–•••

3

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28 29 30 31

94

˜ i }5i =1

Fig. 1. (a), (c) The estimated natural frequencies (vertical axis) {ω

of the Rössler oscillators vs. the natural frequencies of the isolated simulation (horizontal axis).

˜ (θ) estimated by the present method (solid line) is compared with that estimated by the perturbation method (dashed line). All-to-all (b), (d) Interaction function H connections were used to simulate the data in (a) and (b), whereas two defects were considered in the connection matrix to simulate the data in (c) and (d).

32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

95 96 97 98

ω˜ i = 0, a j = b j = 0. To avoid the over-fitting problem, the Fourier

order of D = 1, which represents the simplest form of interaction with smooth oscillator property, was used. This first order approximation has been also justified by the cross-validation test [26]. The convergence property of the multiple-shooting was excellent; few iterations of the Newton procedure gave good estimates. Fig. 1(a), (b) show the estimated natural frequencies of the individual oscillators and the interaction function. 95% confidence intervals indicated by the error-bars were computed from the inverse of the Hessian matrix of the squared error function, based on the assumption that the phase data contain uncorrelated observational noise [27]. The estimated natural frequencies were located on a diagonal line with the original frequencies computed from each of the Rössler oscillators in isolated condition. Moreover, the ˜ (θ) was in a very good agreeestimated interaction function H ment with that estimated by applying the perturbation method to a single isolated Rössler oscillator [2,28]. The estimated connection matrix { T˜ i , j } represented all-to-all connectivity, where all matrix elements were nearly equal to unity (Fig. 2(a)). Next, we introduced two defects, i.e., T 1,2 = T 5,4 = 0, to the connection matrix and estimated the phase dynamics. Although the natural frequencies and the interaction function were estimated quite well in Fig. 1(c), (d), they were slightly deviated from the true ones compared to Fig. 1(a), (b). This is due to the inaccurate assumption of all-to-all connectivity utilized in the first step of the estimation process. The estimation result of the connection matrix, on the other hand, correctly indicated the defects as small values as T˜ 1,2 = 0.163 ± 0.060 and T˜ 5,4 = 0.255 ± 0.0630, whereas other matrix elements pointed almost unity (Fig. 2(b)). Comparable results were obtained in the case that there exist five defects at T 1,2 = T 2,1 = T 3,5 = T 4,3 = T 5,4 = 0, which were detected again as small values as T˜ 1,2 = 0.072 ± 0.076, T˜ 2,1 = 0.0 ± 0.073, T˜ 3,5 = 0.0 ± 0.085, T˜ 4,3 = 0.119 ± 0.077, and T˜ 5,4 = 0.161 ± 0.078

(Fig. 2(c)). Even for a network of N = 7 oscillators having three defects at T 1,2 = T 4,3 = T 6,7 = 0, the corresponding connections were estimates as T˜ 1,2 = 0.125 ± 0.086, T˜ 4,3 = 0.286 ± 0.084, and T˜ 6,7 = 0.270 ± 0.085 (Fig. 2(d)). These results imply that our approach is quite efficient for detecting the defects in the network topology. Accuracy of estimating the connectivity depends upon (1) coupling strength, (2) network size, (3) noise level, (4) defect ratio, (5) length of the measured time series, and (6) threshold to determine the presence of connection. Such dependence was examined in the following. For ouranalysis, the estimation error was evalN N ˆ uated as e = 100 N ( N1−1) i =1 j =1, j =i | T i , j − T i , j | [%], where the

estimated connectivity was digitized as Tˆ i , j = 0 for T˜ i , j < T th and Tˆ i , j = 1 otherwise (T th : threshold). For each setting, 50 instances of connection matrices { T i , j } were randomly generated and the average and the standard deviation of the estimation errors were plotted. Fig. 3(a) shows dependence of the estimation error on the coupling strength C ∈ [0.01, 0.1]. The network size and the number of defects were set to N = 3 and β = 1, where the onset of synchronization was found around C ≈ 0.04. As the coupling strength was increased, the estimation error slowly increased. As discussed in [17,29], this is due to the network dynamics close to the onset of synchronization. When the network is synchronized, phase differences between the oscillators are locked to a constant value θ = const, which provides only poor information about the interaction of the oscillators. For precise estimation of the connectivity, the coupling strength should be small enough to ensure non-synchronized network dynamics. Alternative way of overcoming the problem of synchronized data is to apply external perturbations to destroy the network synchronization and to utilize the relaxation process for the parameter estimation. Next, Fig. 3(b) shows dependence of the estimation error on the network size, which was varied as N = 3, 4, . . . , 10. The number of

99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:PLA

4

AID:21972 /SCO Doctopic: Nonlinear science

[m5Gv1.5; v 1.96; Prn:15/05/2013; 16:08] P.4 (1-6)

I.T. Tokuda et al. / Physics Letters A ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28 29 30 31

94

Fig. 2. Estimation results of the connection matrix { T˜ i , j } ordered in accordance with the connection number i · N + j. The connections corresponding to defects are located on the left side separated from other connections by the dashed line. (a) N = 5. All-to-all connection in original matrix { T i , j }. (b) N = 5. Case of two defects in { T i , j }. (c) N = 5. Case of five defects in { T i , j }. (d) N = 7. Case of three defects in { T i , j }.

32

35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

96 97 98

33 34

95

99

the defects was set to β = 0.2 · N, whereas the coupling strength was increased from C = 0.01 to C = 0.025 so as to locate the data in a non-synchronized regime. These two parameters were adjusted in way that several system sizes can be investigated. As the network size was increased, the estimation error increased. Up to the size of N = 10, the estimation was still reliable with enlarged error bars. Further increase in the network size, however, would eventually lead to unreliable estimates. Dependence of the estimation error on the observational noise level σ ∈ [0, 1] is plotted in Fig. 3(c). Here, independent Gaussian noise N (0, σ ) was added to the phase signals {θi (t )} used for the data fitting. The network size, the defect number, and the coupling strength were set to N = 3, β = 1, and C = 0.005. The estimation error was less than 20% until the noise level of σ = 0.5, which is considered relatively large as a usual condition of experimental measurement. A larger noise, however, increased the estimation error further to 40%, which would eventually reach to a chance level. Although the noise in real environment is usually correlated both in time and between the oscillators, only uncorrelated noise was considered here. Nevertheless, the uncorrelated noise should provide a first step to reveal the main effect of noise on the estimation results. Fig. 3(d) shows dependence of the estimation error on the number of defects β = 1, . . . , N ( N − 1). The network size and the coupling strength were set to N = 5 and C = 0.005. As the number of defects was increased, the estimation error increased due to the violation of all-to-all connectivity assumed in the estimation of the interaction function. The estimation error was, however, less than 20% up to the defect rate of 50%, implying that the present approach is robust enough to be applied to connection matrices with many defects. Fig. 3(e) shows dependence of the estimation error on the length of the time series, i.e., ( M − 1)t. The network size, the

defect number, and the coupling strength were set to N = 4, β = 2, and C = 0.01. As the data length increased, the estimation error decreased. For time series longer than the length of 2000, the estimation error becomes less than 5%. Since natural periods of individual Rössler oscillators are approximately 5.8, 350 cycles of oscillations are included within this data length. Such measurements should be possible for physical or biological systems having a good stationary property. It should also noted that a convergence property of the algorithm was rather independent of the length of the times series. Even for a long time series, few iterations of the Newton procedure converged to a final solution. Finally, Fig. 3(f) shows dependence of the estimation error on the threshold, T th . The network size, the defect number, and the coupling strength were set to N = 4, β = 2, and C = 0.01. For too small or too large setting of the threshold, the estimation error increased. As far as the threshold is chosen as an intermediate value, on the other hand, the estimation error was not too sensitive to a small difference in the threshold. This demonstrates robustness of the present results against the setting of the threshold.

100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119

4. Application to electrochemical oscillators

120 121

We tested the applicability of our technique to experimental data with a small network of electrochemical oscillators. As a preliminary test, our method was first applied to simulated data from a mathematical model of the electrochemical oscillators [30] and we have obtained essentially the same results with those of the Rössler oscillators. Then, the experiments were carried in a standard electrochemical cell containing 3 mol/dm3 sulfuric acid at 10 ◦ C. Three 1 mm diameter Ni electrodes serve as the working, a Hg/Hg2 SO4 /sat. K2 SO4 as reference, and a Pt rod as the counter electrode. The experimental setup is identical to that applied in our previous studies [13,23]. The three Ni electrodes are held at

122 123 124 125 126 127 128 129 130 131 132

JID:PLA AID:21972 /SCO Doctopic: Nonlinear science

[m5Gv1.5; v 1.96; Prn:15/05/2013; 16:08] P.5 (1-6)

I.T. Tokuda et al. / Physics Letters A ••• (••••) •••–•••

5

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43 44 45 46

Fig. 3. (a) Dependence of the estimation error on the coupling strength C ∈ [0.01, 0.1] used to generate the simulation data. For 50 random realizations of the connection matrix, the average and the standard deviations of the estimation errors are plotted. (b) Dependence of the estimation error on the network size N = 3, 4, . . . , 10. (c) Dependence of the estimation error on the observational noise σ ∈ [0, 1]. (d) Dependence of the estimation error on the defect rate β/ N ( N − 1). (e) Dependence of the estimation error on the length of the measured time series ( M − 1)t. (f) Dependence of the estimation error on the threshold, T th , which judges the presence of connection.

47

50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

110 111 112 113

48 49

109

114

the same constant potential of 1080 mV by a potentiostat. The currents of the electrodes, proportional to the oscillating rate of electrodissolution, are measured by digitizing the potential drop on individual resistors (R 1 = 950 , R 2 = 1000 , and R 3 = 1050 ) attached to the electrodes. (We applied slightly different individual resistors to induce reproducibly different natural frequencies of the oscillators.) Coupling between the electrodes is induced by cross resistances (R 1,2 and R 2,3 ) between the nickel wires. The oscillations of current occur through a Hopf bifurcation at a potential of about 1060 mV. The oscillators at the operating conditions are considered smooth with first harmonic components dominating the phase interaction function [31]. From the network of N = 3 electrochemical oscillators, multivariate current data { y i (t )}3i =1 were measured and the phases were reconstructed in a similar manner as for the simulation data. Since first harmonic components dominate the interaction function of the experimental system, the order of D = 1 was used. Two instances of network topology were considered. The first one

has two defects at T 1,3 = T 3,1 = 0, whereas the second one has four defects at T 1,3 = T 2,3 = T 3,1 = T 3,2 = 0. The coupling strength was tuned in such a way that the electrochemical system generates non-synchronized dynamics with phase slipping behavior (see Fig. 4). The connection matrix and the natural frequencies were estimated for the two instances. As shown in Fig. 4, the estimated natural frequencies are within 1 mHz of the actual natural frequencies, which is about 6% of the typical natural frequency differences (16 mHz). This 6% deviation is sufficiently accurate for reconstruction of synchronization patterns in the experimental setting. The connection matrix was estimated as follows:

⎛

T˜ 1,2

⎝ T˜ 2,1 T˜ 3,1 T˜ 3,2  =

⎞

116 117 118 119 120 121 122 123 124 125 126

T˜ 1,3 T˜ 2,3 ⎠

1.0 ± 0.08 0.27 ± 0.10 1.0 ± 0.11 1.0 ± 0.09 0.18 ± 0.10 1.0 ± 0.10

115

127 128 129



130

,

131 132

JID:PLA

AID:21972 /SCO Doctopic: Nonlinear science

[m5Gv1.5; v 1.96; Prn:15/05/2013; 16:08] P.6 (1-6)

I.T. Tokuda et al. / Physics Letters A ••• (••••) •••–•••

6

One of the severe limitations of the present parameter estimation may arise from the system size. As the system size increases, the size of the connection matrix grows as ∼ N ( N − 1). According to our study, reliable results have been obtained for a system size of up to N = 10. The system size could be too small for single-cell network problems in neuroscience, where often network sizes consisting of hundreds of neurons arise. However, the method could be used for identification of connectivity between relatively small number of brain regions or in network motifs of brain circuitries [34]. In our future work, comparative study with other methodologies will be carried out to clarify advantages as well as disadvantages of our approach. One of our future challenges is to apply the present technique to a network of biological oscillators such as heart-respiratory interaction [35] or fetal-maternal heart interaction [36]. The phase models built from such experimental data would facilitate an essential description and new understanding of the biological functions.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

Acknowledgements

21

This work was partially supported by Grant-in-Aid for Scientific Research (C) (No. 23560446) from Japan Society for the Promotion of Science (JSPS) and by the National Science Foundation under Grant No. CHE-0955555.

23 24 25

28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

Fig. 4. Experiments: Reconstruction of network topology of electrochemical oscillator networks. Left column: Three locally coupled oscillators with two defects. R 1,2 = 37 k, R 2,3 = 42 k. Right column: Three oscillators with four defects. R 1,2 = 75 k. Top row: Phase difference between the elements in the network. Middle row: Estimated vs. natural frequencies. Bottom row: Reconstructed connectivity matrix.

⎛

T˜ 1,2

⎝ T˜ 2,1 T˜ 3,1 T˜ 3,2  =

⎞

T˜ 1,3 T˜ 2,3 ⎠ 0.77 ± 0.19

1.0 ± 0.28 0.05 ± 0.23 0.19 ± 0.25

0.0 ± 0.24 0.13 ± 0.19

.

We see that the defects were precisely identified as small values in the estimated connection matrix, whereas the other matrix elements pointed almost unity. 5. Conclusions and discussions

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 86 87

22

27

68

85

20

26

67

To conclude, a simple yet efficient methodology has been presented for estimating the connection matrix of oscillator networks from multivariate time series. Our methodology was efficient for the data simulated from a network of up to 10 limit cycle oscillators. To ensure reliable estimates, the network dynamics should be in a non-synchronized state. Robustness of the estimation against observational noise has been also confirmed. The practical applicability to experimental systems was demonstrated by the accurate detection of the defects in the measurement data from the electrochemical oscillators. It should be noted that the present approach provides estimates not only for the connection matrix but also for the natural frequencies and the interaction function of the phase dynamics. This provides a valuable insight into the detailed dynamics of the observed oscillator networks. Because of the all-to-all connectivity assumed for the estimation of the natural frequencies and the interaction function, the accuracy can be lowered for a network of sparsely connected oscillators. We may, however, expect that in some systems the interaction function as well as the natural frequencies can be measured a priori, e.g., by conventional techniques [28,31–33]. Use of such prior estimates would largely reduce the difficulty of our technique and enhance its wider applications to networks having many defects.

88 89 90 91 92

References

93 94

[1] A.T. Winfree, The Geometry of Biological Time, Springer, New York, 1980. [2] Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer, Berlin, 1984. [3] A. Pikovsky, M. Rosenblum, J. Kurths, Synchronization – A Universal Concept in Nonlinear Sciences, Cambridge University Press, Cambridge, 2001. [4] D. Watts, S. Strogatz, Nature 393 (1998) 440. [5] S. Yamaguchi, H. Isejima, T. Matsuo, R. Okura, K. Yagita, M. Kobayashi, H. Okamura, Science 302 (2003) 1408. [6] S.J. Aton, E.D. Herzog, Neuron 48 (2005) 531. [7] T. Schreiber, Phys. Rev. Lett. 85 (2000) 461. [8] M.G. Rosenblum, A.S. Pikovsky, Phys. Rev. E 64 (2001) 045202. [9] M.C. Romano, M. Thiel, J. Kurth, C. Grebogi, Phys. Rev. E 76 (2007) 036211. [10] S. Hempel, A. Koseska, J. Kurths, Z. Nikoloski, Phys. Rev. Lett. 107 (2011) 054101. [11] B. Schelter, M. Winterhalder, R. Dahlhaus, J. Kurths, J. Timmer, Phys. Rev. Lett. 96 (2006) 208103. [12] J. Nawrath, M.C. Romano, M. Thiel, I.Z. Kiss, M. Wickramasinghe, J. Timmer, J. Kurths, B. Schelter, Phys. Rev. Lett. 104 (2010) 038701. [13] M. Wickramasinghe, I.Z. Kiss, Phys. Rev. E 83 (2011) 016210. [14] J. Runge, J. Heitzig, V. Petoukhov, J. Kurths, Phys. Rev. Lett. 108 (2012) 258701. [15] M. Timme, Phys. Rev. Lett. 98 (2007) 224101. [16] Z. Levnajic, A. Pikovsky, Phys. Rev. Lett. 107 (2011) 034101. [17] I.T. Tokuda, S. Jain, I.Z. Kiss, J.L. Hudson, Phys. Rev. Lett. 99 (2007) 064101. [18] B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, Phys. Rev. E 77 (2008) 066205. [19] B. Kralemann, A. Pikovsky, M. Rosenblum, Chaos 21 (2011) 025104. [20] K.A. Blaha, A. Pikovsky, M. Rosenblum, M.T. Clark, C.G. Rusin, J.L. Hudson, Phys. Rev. E 84 (2011) 046201. [21] K. Wiesenfeld, J.W. Swift, Phys. Rev. E 51 (1995) 1020. [22] K. Wiesenfeld, C. Bracikowski, G. James, R. Roy, Phys. Rev. Lett. 65 (1990) 1749. [23] I.Z. Kiss, Y.M. Zhai, J.L. Hudson, Science 296 (2002) 1676. [24] D. Hansel, H. Sompolinsky, Phys. Rev. Lett. 68 (1992) 718. [25] E. Baake, M. Baake, H.G. Bock, K.M. Briggs, Phys. Rev. A 45 (1992) 5524. [26] M. Stone, J. Royal Stat. Soc. B 36 (1974) 111. [27] W. Horbelt, J. Timmer, M.J. Bünner, R. Meucci, M. Ciofini, Phys. Rev. E 64 (2001) 016222. [28] H. Sakaguchi, S. Shinomoto, Y. Kuramoto, Prog. Theor. Phys. 77 (1987) 1005. [29] I. Tokuda, J. Kurths, E. Rosa Jr., Phys. Rev. Lett. 88 (2002) 014101. [30] Y. Zhai, I.Z. Kiss, J.L. Hudson, Ind. Eng. Chem. Res. 43 (2004) 315. [31] I.Z. Kiss, Y.M. Zhai, J.L. Hudson, Phys. Rev. Lett. 94 (2005) 248301. [32] R.F. Galan, G.B. Ermentrout, N.N. Urban, Phys. Rev. Lett. 94 (2005) 158101. [33] J. Miyazaki, S. Kinoshita, Phys. Rev. Lett. 96 (2006) 194101. [34] O. Sporns, R. Koetter, PLoS Biol. 2 (2004) e369. [35] C. Schäfer, M.G. Rosenblum, H.-H. Abel, J. Kurths, Phys. Rev. E 60 (1999) 857. [36] P. Van Leeuwen, D. Geue, M. Thiel, D. Cysarz, S. Lange, M.C. Romano, N. Wessel, J. Kurths, D.H. Grönemeyer, Proc. Natl. Acad. Sci. USA 106 (2009) 13661.

95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132