Dynamical Bayesian Inference of Time-evolving Interactions: From a Pair of Coupled Oscillators to Networks of Oscillators

1.1 MB·English

by  Andrea Duggento

DynamicalBayesianInference ofTime-evolvingInteractions: From aPairofCoupled Oscillatorsto
Networks ofOscillators
AndreaDuggento1, Tomislav Stankovski2, Peter V. E. McClintock2, and Aneta Stefanovska2∗
1 Medical PhysicsSection, Faculty of Medicine, Tor Vergata University, Rome, Italyand
2 Department of Physics, Lancaster University, Lancaster, LA1 4YB, United Kingdom
(Dated:September24,2012)
Livingsystemshavetime-evolvinginteractionsthat, untilrecently, couldnot beidentifiedaccuratelyfrom
recordedtimeseriesinthepresenceofnoise. Stankovskietal.(Phys.Rev.Lett. 109024101;2012)introduced
a method based on dynamical Bayesian inference that facilitates the simultaneous detection of time-varying
2
synchronization,directionalityofinfluence,andcouplingfunctions.Itcandistinguishunsynchronizeddynamics
1
fromnoise-inducedphaseslips. Themethodisbasedonphasedynamics,withBayesianinferenceofthetime-
0
evolvingparametersbeingachievedbyshapingthepriordensitiestoincorporateknowledgeofprevioussamples.
2
Wenowpresentthemethodindetailusingnumerically-generateddata,datafromananalogelectroniccircuit,

p andcardio-respiratorydata.Wealsogeneralizethemethodtoencompassnetworksofinteractingoscillatorsand
e thusdemonstrateitsapplicabilitytosmall-scalenetworks.
S

1 PACSnumbers: 02.50.Tt,05.45.Xt,05.45.Tp,87.10.-e,87.19.Hh
2

] I. INTRODUCTION directionality[21–23]. Thesetechniquesprovidemeasuresof
n theamountofinformationinameasuredsignal,orthecausal
a
relationshipsbetweenmeasuredsignalsand,indoingso,they
- Systemsofinteractingoscillatorsareubiquitousinscience.
a infereffectratherthanmechanism.
In the commoncase where the naturalfrequenciesor ampli-
t
a tudesorinter-oscillatorcouplingsaretime-varying,theypose Inthispaperwedescribeindetailandfurtherextendare-
d acontinuingchallengetothetime-seriesanalystwhoendeav- centlyintroducedmethod[24] basedondynamicalBayesian
s. ours to understand the underlying system from the signal(s) inference. As we demonstratebelow, it enjoysmany advan-
c it creates. Oversimplifications of hypotheses are often used tagesoverearlierapproaches. Oneoftheseisthatitdoesnot
i
s to render the problem more tractable, but can all too easily requiretheobservabletofillthedomainoftheprobabilityden-
y result in a failure to describe phenomena that are in fact of sity functionatequilibrium: itcanthusprovidethesame in-
h
central importance – given that the strength, direction and formationfromasmallfractionofthedatathatisrequiredby
p
[ functionalrelationships that define the nature of the interac- the transfer entropy or Granger causality approaches. Addi-
tions can cause qualitatively new states to appear or disap- tionally,thedynamicalapproachhastheadvantageofcharac-
1 pear. Time-variabilityof this kind is especially importantin terizingthesystemcompletely(notonlyintermsofinforma-
v biologicalapplications,thoughitisbynomeansrestrictedto tion measures). Thus, from the inferred dynamics one can
4
biology. deduce self-consistently any information that is of interest,
8
be it coupling functions, or synchronization, or causality, or
6 Intheabsenceoftime-variability,therearemanydifferent
equilibriumdensities,includingtheequationsofmotion. We
4 methodsavailable[1–4]fordetectingandquantifyingthecou-
9. plingsanddirectionality(dominantdirectionofinfluence)be- discussindetailthetheoreticalbackground,thetechnicalas-
pects and limitations of the algorithms, and we demonstrate
0 tween oscillators based, especially, on the analysis of phase
thewideapplicabilityofthemethodbyconsiderationofsev-
2 dynamics. Approaches to the detection of synchronization
1 have mostly been based on the statistical properties of the eralexamples.
v: phasedifference[5–8]. Theinferenceofanunderlyingphase Thecouplingfunctionsareofparticularimportance. Their
i modelhasbeen used as the functionalbasis for a numberof form is uniquely able to describe the functional laws of in-
X techniquestoinferthenatureofthephase-resettingcurves,in- teraction between the oscillators. Earlier theoretical studies
r teractionsandstructuresofnetworks[9–14]. However,these have included the work of Kuramoto[25] and Winfree [26],
a
techniques inferred neither the noise dynamics nor the pa- whichused a functiondefinedeither by the phase difference
rameterscharacterizingthenoise. Anadditionalchallengeto or by both phases, and of Daido [27, 28] and Crawford[29]
these methods can be the time-varying dynamics mentioned who used a more general form in which the coupling func-
above. Inaseparatelineofdevelopment,Bayesianinference tion was expanded in Fourier series. Other methods for in-
wasappliedtoanalysethesystemdynamics[15–20],thereby ference of the coupling functions have also been suggested
opening the door to inference of noisy time-evolving phase [10,11,13,30]. Thetechniquedescribedbelowgoesbeyond
dynamics. Methods based on transfer entropy and Granger allofthesebecauseitisabletofollowthetime-variabilityof
causalityhaveageneralitythathasfacilitatedanumberofap- the coupling functionsand hence can revealtheir dynamical
plications, including inference of the coupling strength and characterwhereitexists.
We will also show how the technique can readily be
extended to encompass networks of interacting oscillators.
Theseformalargeandimportantgroupofphysicalsystems,
∗Electronicaddress:[email protected] including neural networks [9, 31, 32], electrochemical sys-
2
tems [10, 33], crowd synchrony on the Millennium bridge, probabilitydistribution(seeSec. IIBforstationarydy-
andnetworksoffireflies[34]. Thelargescaleofthenetworks namics,andSec. IICfortime-varyingdynamics).
canintroduceahighercomplexity,bothinstructureandfunc-
3. Integrationoftheprobabilitythatthisparametersetlies
tional behavior. For example in neuronal networks, the ex-
insidetheArnoldtonguedefiningsynchronization.This
istenceofspatialandspatial-temporalcorrelations,collective
effectivelyyieldsthecumulativeprobabilityofthesyn-
or partially collective (clustering) behavior, synchronization
chronizationstateofthedynamics(seeSec. IIIA).
or desynchronization, and time-variability has been reported
[31, 32, 35, 36]. In such cases, and given the kinds of phe-
4. Useoftheparameterinformationasobtainedinstep2
nomenontobestudied,thereisanincreasingneedforpower-
tocreateadescriptionoftheinteractions,leadingtode-
fultechniquesthatcaninferthetime-varyingdynamicsofthe
tection of the predominantdirectionality and coupling
oscillatorynetworks.
functionestimationamongtheoscillators(seeSec.IV).
In Sec. II we provide details about the phase decomposi-
tions,theimplementationoftheBayesianframeworkandhow
thetime-varyinginformationispropagated.Thesynchroniza-
A. TruncatedFourierSeries
tion detectionthrougha maprepresentationof the phase dy-
namicsisdiscussedinSec.III,whilethemethodfordescrib-
The periodic behaviour of the system suggests that it can
ing the interactions is demonstrated in Sec. IV. Before the
appropriatelybe described by a Fourier decomposition. De-
methodisapplied,weconsiderinSec.Vsomeimportanttech-
composingbothf andg inthiswayleadstotheinfinitesums
nical aspects and limitations. The wide applicability of the i i
method is demonstrated in Sec. VI, through the analysis of ∞
time-seriesfromnumericalphase and limit-cycleoscillators, f (φ ) = c˜ sin(kφ )+c˜ cos(kφ )
i i i,2k i i,2k+1 i
analoguesimulationandcardio-respiratoryinteractions. The k=X−∞
generalization of the approach to networks of oscillators, as and
exemplifiedbytwo numericalexamples,is presentedin Sec. ∞ ∞
VII.Finally,wesummariseanddrawconclusionsinSec.VIII. gi(φi,φj) = c˜i;r,sei2πrφiei2πsφj. (2)
Thealgorithmusedforthedetectionofsynchronizationisde- s=X−∞r=X−∞
scribedinAppendixA.
It is reasonable to assume that, in most cases, the dynamics
willbewell-describedbyafinitenumberK ofFourierterms,
sothatwecanrewritethephasedynamicsofEq.(1)asafinite
II. PHASE-DYNAMICSDECOMPOSITION
sumofbasefunctions
Consider an N-dimensional oscillator dx/dt = f(x(t)) K
whosesolutionf admitsalimitcycle. Suchanoscillatorcan φ˙ = c(l)Φ (φ ,φ )+ξ (t), (3)
l k l,k 1 2 l
usuallyberepresentedbyaconstantphasevelocityφ˙ =ωand k=X−K
avectorcoordinatethatdefinesthelimitcycleasafunctionof
thephaseφ: r r(φ). wherel=1,2,Φ = Φ =1,c(l) = ω ,andtherestofΦ
≡ 1,0 2,0 0 l l,k
When two such oscillators mutually interact sufficiently andc(l) aretheK mostimportantFouriercomponents.
k
weakly,theirmotioniscommonlyapproximatedjustbytheir
phase dynamics [25, 37]. We note that, in general, if we
describe the phase of a system through a generic monotonic B. BayesianInference
changeofvariables,thanthe dynamicalprocesscan bewrit-
tenas
In orderto reconstructthe parametersof Eq. (3) we make
φ˙ =ω +f (φ )+g (φ ,φ )+ξ . (1) extensive use of the approach to dynamical inference pre-
i i i i i i j i
sentedin[18,19]. Inthissectionwebrieflyoutlinethetech-
Eq.(1)explicitlyincludesanoisetermξitoenableittorepre- niqueasadaptedto thepresentcase. Thefundamentalprob-
sentaprocessinarealsystem. Thenoisecanbee.g.awhite lem in dynamical inference can be defined as follows. A 2-
Gaussiannoise ξi(t)ξj(τ) = δ(t τ)Eij,whichmaypos- dimensional(ingeneralL-dimensional)time-seriesofobser-
h i −
siblybespatiallycorrelated. vationaldata = φ φ (t ) (t = nh) is provided,
l,n l n n
The phase-dynamics decomposition technique is highly andtheunknoXwnm{odelpa≡rameters} = c(l),E , areto
modular from the algorithmic point of view, and each mod- M { k ij }
beinferred.
ule will be explained separately in the sections that follow.
Bayesianstatisticsemploysagivenpriordensityp ( )
prior
Theoverallprocedurecomprisesthefollowingsteps: M
that encloses expert knowledge of the unknown parameters,
1. Assumption that the dynamics can be precisely de- together with a likelihood function ℓ( ), the probability
X|M
scribed by a finite number of Fourier terms (see Sec. density to observe φl,n(t) given the choice of the dy-
{ } M
namicalmodel.Bayes’theorem
IIA).
2. Inference,giventhedata,oftheFourierterms,thenoise ℓ( )pprior( )
pX( )= X|M M (4)
amplitude,andtheircorrelationinformofaparameter M|X ℓ( )p ( )d
prior
X|M M M
R
3
then enables calculation of the so-called posterior density Σ , the stationarypointofS canbecalculatedrecursively
prior
pX( )oftheunknownparameters conditionedonthe from
M|X M
observations.
h
miFdo-proiinndteappepnrdoexnimtwathioitnewGhaeursesian noise sources, and in the Eij = N (cid:16)φ˙i,n−ck(i)Φi,k(φ∗·,n)(cid:17)(cid:16)φ˙j,n−c(kj)Φj,k(φ∗·,n)(cid:17),
c(i) =(Ξ−1)(i,l)r(l),
k kw w
φ φ (φ +φ )
φ˙l,n = l,n+1h− l,n and φ∗l,n = l,n 2 l,n+1 rw(l) =(Ξ−pri1or)k(iw,l)c(wl)+hΦi,k(φ∗·,n)(E−1)ijφ˙j,n+
the likelihood is given by a product over n of the probabil- h∂Φl,k(φ·,n),
ityofobservingφl,n+1 ateachtime. Thelikelihoodfunction − 2 ∂φl
is constructed by evaluation of the stochastic integral of the Ξ(i,j) =Ξ (i,j)+hΦ (φ∗ )(E−1) Φ (φ∗ ),
kw priorkw i,k ·,n ij j,w ·,n
noisetermovertime,as
(7)
ξ(1)(t ) ti+1ξ (t)dt=√hHz , (5) where the covariance is Σ=Ξ−1, summation over n from
l i ≡Zti l l 1 to N is assumed and summation over repeated indices
k,l,i,j,wisagainimplicit.
whereH isthe Choleskydecompositionofthepositivedefi-
We note that a noninformative“flat” prior can be used as
nitematrixE, andz isaspatiallyuncorrelatedstandardized
l the limitof an infinitelylargenormaldistribution,by setting
Gaussianvariable. Thejointprobabilitydensityofzl isused Ξ =0andc¯ =0.
tofindthejointprobabilitydensityoftheprocessinrespectof prior prior
Themultivariateprobability X(c,c¯,Ξ)giventhereadout
(φ (t ) φ(t ))byimposingP(φ (t )=det(Jφ)P(ξi), N |
l i+1 − l i l i+1 ξ timeseries = φl,n φl(tn) explicitlydefinestheprob-
where Jφ is the Jacobian term of the transformationof vari- abilitydensXityof{eachp≡arameter}setofthedynamicalsystem.
ξ
ablesthatcanbecalculatedfromEq.(2). Ifthesamplingfre- Because each of them can be discriminated, as belongingor
quencyishighenough,thetimestephtendstozero,andthe notbelongingtothe Arnoldtongueregionwe candefinethe
determinantoftheJacobianJφ canbewell-approximatedby binarypropertys(c(l))= 1,0 ,andcanobtaintheposterior
ξ k { }
theproductofitsdiagonalterms probabilityofthesystembeingsynchronizedornotbyevalu-
atingtheprobabilityofs
det(Jφk(tn)) ∂Φl,k(φ·,n).
ξk(tn) ≈Yl ∂φl psync ≡pX(s=1)=Z s(c)NX(c|c¯,Ξ)dc. (8)
This transformationleads to an extra term in a least squares
Thecomputationofp willbediscussedinSec. III.
likelihood, and the minus log-likelihood function S = sync
lnℓ( )canthusbewrittenas
− X|M
C. Time-varyinginformationpropagation
N−1
S = N ln E + h c(l)∂Φl,k(φ·,n)+
2 | | 2 nX=0(cid:16) k ∂φl ThemultivariateprobabilitydescribedbyNX(c,Σ)forthe
given time series = φ φ(t ) explicitly defines the
+[φ˙i,n−ck(i)Φi,k(φ∗·,n)](E−1)ij[φ˙j,n−c(kj)Φj,k(φ∗·,n)](cid:17), probability densityXof e{achnp≡aramente}r set of the dynamical
(6) system. When the sequential data come from a stream of
measurementsprovidingmultipleblocksof information,one
wheresummationovertherepeatedindicesk,l,i,jisimplicit. applies (7) to each block. Within the Bayesian theorem, the
The log-likelihood (6) is a quadratic form of the Fourier evaluationof the currentdistributionrelies on the evaluation
coefficientsofthephases. Henceifamultivariatepriorprob- ofthepreviousblockofdata,i.e.thecurrentpriordependson
ability is assumed, the posterior probabilityis a multivariate thepreviousposterior. Thustheinferencedefinedinthisway
normaldistributionaswell. is not a simple windowing, but each stationary posterior de-
Thisishighlydesirablefortworeasons:(i)aGaussianpos- pendsonthe historyof theevaluationsfrompreviousblocks
terior is computationally convenient because it guarantees a ofdata.
unique maximum, with the mean vector and covariancema- In classical Bayesian inference, if the system is known to
trix completely characterizing the distribution and giving us benon-time-varying,thentheposteriordensityofeachblock
themostsignificantinformation;(ii)allthemultivariatenor- istakenasthepriorofthenextone: Σn+1 = Σn . Thisfull
prior post
malposteriorscanbeusedagainaspriorsinthepresenceofa propagationof the covariancematrix allowsgoodseparation
newblockofdata, andknowledgeaboutthesystemcaneas- of the noise, and the uncertaintiesin the parameterssteadily
ilybeupdated. Thislastfeatureisessentialforanyreal-time decreasewithtimeasmoredataareincluded.
applicationbecauseitensuresthatthecomplexityofthealgo- Iftime-variabilityexists,however,thispropagationwillact
rithmdoesnotchangewiththelengthoftheinputdata-stream. as a strong constraint on the inference, which will then fail
From[18],andassumingamultivariatenormaldistribution to follow the variations of the parameters. This situation is
asthepriorforparametersc(l),withmeansc¯,andcovariances illustrated in Fig. 1(a) [53]. In such cases, one can consider
k
4
the processes between each block of data to be independent
(i.e.Markovian).Therecannotthenbeanyinformationprop- 4
agationbetweentheblocksofdata,andeachinferencestarts 3.5
fromthe flat distributionΣnpr+io1r = 0. The inferencecan thus ε(t)2 3
followmorecloselythetime-variabilityoftheparameters,but
theeffectofnoiseandtheuncertaintyoftheinferencewillof 2.5 (a) (b) (c)
coursebemuchlarger,asshowninFig.1(b). 0 400 8000 400 8000 400 800
Time [s] Time [s] Time [s]
Where the system’s parameters are time-dependent, we
may assume that their probability diffuses normally accord-
Figure 1: Inference of a rapidly time-varying coupling parameter
ingly to the known diffusion matrix Σ . Thus, the proba-
diff
fromcouplednoisyoscillators(12). Thegraylinesrepresenttheac-
bilitydensityofthe parametersistheconvolutionofthe two
tualparameterinthenumericalsimulation,whereastheblacklines
normalmultivariatedistributions,Σ andΣ
post diff indicatethetime-varyingparameterinferredfromtheresultanttime
Σnpr+io1r =Σnpost+Σndiff. stieornie,sΣ, nfo+r1: =(a)∞fu;llapnrdop(ca)gaptrioopna,gΣatnpiro+ionr1fo=r tiΣmnpeos-tv;a(rby)inngoprporcoepsasgeas-,
ThecovariancematrixΣdiff expressesourbeliefaboutwhich Σn+1 =prioΣrn +Σn .
part of the dynamical fields that define the oscillators has prior post diff
changed, and the extent of that change. Its elements are
(Σ ) =ρ σ σ ,whereσ isthestandarddeviationofthe
diff i,j ij i j i
diffusionoftheparameterc afterthetimewindowt thathas
i w III. SYNCHRONIZATIONDETECTION
elapsed from the first block of information to the following
one. ρ isthecorrelationbetweenthechangeoftheparame-
ij
Itisimportanttonotethatfinitenoisecaninducephaseslips
tersc andc (withρ =1). Inrelationtothelatter,aspecial
i j ii
inasystemthatwouldbesynchronizedinthenoiselesslimit.
exampleof Σ will be considered: we assume thatthere is
diff
Ratherthanfocusingonthe presenceandstatistics ofphase-
nocorrelationbetweenparameters,i.e.ρ =0,andthateach
ij
slips, we propose to detect synchronization from the nature
standarddeviationσ isaknownfractionoftheparameterc :
i i
ofthephase-slipitself. Anovelfeatureofthepresentstudyis
σ = p c (wherep indicatesthatpisreferredtoawindow
i w i w
thatitproposesevaluationoftheprobabilitythattheequations
oflengtht ). Itisimportanttonotethatthisparticularexam-
w
drivingthe dynamicsare intrinsically synchronizedand thus
pleisactuallyrathergeneralbecauseitassumesthatallofthe
parameters(fromtheΣn diagonal)canbeofatime-varying ofwhetheranyphase-slipsthatmaypossiblybeobservedare
post
dynamics-relatedornoise-induced.
nature – which corresponds to the inference of real (experi-
mental)systemswithaprioriunknowntime-variability. After performing the inference, one can use the recon-
There are two obvious limits in modeling the knowledge structed parameters, derived in the form of a multivariate
assumed with respect to possible time variation of param- normal distribution X(c,Σ), to study the interactions be-
N
eters. The first of these is to assume no time-variability: tween the oscillators under study. In general, the border of
in this case the full information propagation matrix is used, theArnoldtonguemaynothaveananalyticformand,evenif
n+1 = n . If the assumption proves wrong, the in- ithas,theintegralmaynothaveananalyticsolutionbutmust
prior post be evaluated numerically. In practice, we estimate p nu-
Pferred paramPeters may accumulate a bias when the real sys- sync
merically, sampling from the parameter space many realiza-
tem varies in time. The other limit is to assume each time
windowto be completelyindependentof the previoussignal tions {c(kl)}m, where m labels each parameter vector tested.
history. In this case no propagation is used, n+1 = , Foreverysetofcwecomputes(cm)numerically. Letusas-
prior ∞ sume for now that s(c ) is given. To find p with arbi-
(i.e.Ξn+1 =0),andthereisnobiasbut,becauPsemuchinfor- m sync
prior trary precision, it is enough to generate a number M of pa-
mationisforgotten,theprobabilityoftheinferredparameters
rameters c = c(l) , with m = 1,...,M sampled from
has a large covariance matrix. An optimal assumption must m { k }m
lies in between these two limits: npr+io1r = npost+ ndiff. NX(c|c¯,Ξ),sincepsync =limM→∞ M1 Mms(cm).
If a diffusionmatrix is assumed, wPe allow thPe methodPsome However, this 2K-dimensional intPegration quickly be-
freedomforthetime-variabilitytobefollowed,whilerestrict- comesinefficientwithanincreasingnumberofFouriercom-
ing it to be unbiased. The amount of variability is part of ponents. Moreover,aswewilldiscussinSec.IIIA,thecom-
themodel,likethenumberoffreeparametersinanystandard putationtimeofthevariables(cm)isnotinsignificant.Onthe
method. Fig.1illustratesthetwoextremelimits,andapossi- other hand, if the posterior probabilitypX is sharply peaked
bletrade-off. TheinferenceinFig.1(c)demonstratesthatthe aroundthe meanvaluec¯, thenpsync willbe indistinguishable
time-variabilityiscapturedcorrectlyandthattheuncertainty froms(c¯),andtheevaluationofs(c¯)willsuffice.
isreducedbecausemoredatahavebeenincluded.
If one knows beforehandthat only one parameter is vary-
ing(or,atmost,asmallnumberofparameters),thenΣ can A. SynchronizationDiscriminationandmaprepresentation
diff
becustomizedtoallowtrackingofthetime-variabilityspecif-
ically of that parameter. This selective propagation can be We now illustrate a simple numerical technique to recog-
achievedif,forexample,notallbutonlytheselectedcorrela- nize whether a coupled phase oscillator system is synchro-
tionρ fromthediagonalhasanon-zerovalue. nized,ornot. Thetechniqueitselfamountstoasimplecheck
ii
5
functions f (φ ) and g (φ ,φ ) that are linked to influences
i j i i j
betweentheoscillators.
Onecan seek to determinethe propertiesthatcharacterize
the interaction in terms of a strength of coupling, predomi-
nantdirectionofcoupling,orevenbyinferenceofacoupling
function. Theanalysisofinformationpropagationallowsin-
=0
ference of the time-varying dynamics, and the interactions’
propertiescanbetracedintimeaswell. Thisisespeciallyim-
Figure2: Torusrepresentationofthephasedynamics,withtoroidal portant for the inference of open interacting oscillatory pro-
coordinate ζ(φ1(t),φ2(t)) and polar coordinate ψ(φ1(t),φ2(t)). cesses wherethe time-variabilityof the interactionscan lead
ThewhitecircledenotesthePoincarécrosssection.
to transitions between qualitatively different states, such as
synchronization[37]oroscillatingdeath[38].
The coupling amplitude quantifies the total influence be-
tweentheoscillatorsinaparticulardirection,e.g.howmuch
by numerical integration of the system of ordinary differen-
the dynamicsof the first oscillator affects the dynamicalbe-
tial equationdefinedbyEq. (1) throughonecycleof the dy-
havior of the second oscillator (1 2). Depending on
namics,andtestingwhetherthe1:1synchronizationcondition →
whether the coupling is in only one direction, or in both di-
ψ(t) = φ (t) φ (t) <K isalwaysobeyed.
| Let| us|as1sum−e w2e ar|e observing motion on the torus T2 rections,wespeakofunidirectionalorbidirectionalcoupling,
respectively.Intheinferentialframeworkthatwepropose,the
definedbythetoroidalcoordinateζ(φ (t),φ (t))=(φ (t)+
1 2 1 coupling amplitudes are evaluated as normalized measures,
φ (t))/2,andthepolarcoordinateψ(t).
2 basedontheinteractingparametersinferredfromthecoupling
For assessment of possible 1:1 synchronization the phase
base functionsq (φ ,φ ). The quantificationiscalculatedas
differenceψ(t)willbedefinedasψ(φ (t),φ (t)) = φ (t) i i j
1 2 1 − aEuclidiannorm:
φ (t). Fig.2providesaschematicrepresentationofthephase
2
dynamicsonthetorus. LetusconsideraPoincarésectionde-
ǫ = q (φ ,φ ) c2+c2+...
fined byζ = 0 andassume thatζ(t)/dt|ζ=0 > 0 forany ψ. 21 k 1 1 2 k≡q 1 3 (9)
Thismeansthatthedirectionofmotionalongthetoroidalco-
ǫ = q (φ ,φ ) c2+c2+...,
ordinateisthesameforeverypointofthesection. Ideallywe 12 k 2 1 2 k≡q 2 4
wouldfollowthe time-evolutionofeverypointandestablish
where the odd inferred parameters are assigned to the base
whetherornotthereisaperiodicorbit;ifthereisone,andif
functionsq (φ ,φ )forthecouplingthatthesecondoscillator
itswindingnumberiszero, thenthesystem issynchronized. 1 1 2
imposesonthefirst(ǫ : 2 1),andviceversa(ǫ : 1
Ifsucha periodicorbitexists, thenthereisatleastoneother 21 → 12 →
2).
periodic orbit, with one of them being stable and the other
Thedirectionalityofcoupling[2]oftenprovidesusefulin-
unstable.
formationabouttheinteractions.Itisdefinedasnormalization
Thesolutionofthedynamicalsystemoverthetorusyields
aboutthepredominantcouplingamplitude
amapM : [0,2π] [0,2π]thatdefines,foreachψ onthe
n
→
tPoorionicdaarlécsoeocrtdioinna,tteheψnext=phaMse(ψψn+).1Faifgte.r3o(bn)e,(cci)rciulliutsotrfatthees D = ǫ12−ǫ21. (10)
n+1 n ǫ +ǫ
12 21
the map M as evaluated computationally in two situations,
correspondingto nosynchronization,or synchronization,re- IfD (0,1]thefirstoscillatordrivesthesecond(1 2),or
∈ →
spectively. ifD [ 1,0)thesecond(2 1)drivesthefirst. Thequan-
∈ − →
The map M is continuous, periodic, and has two fixed tified valuesofthe couplingstrengthsǫ orthedirectionality
i
points (one stable and one unstable) if and only if there is a D representmeasuresofthecombinedrelationshipsbetween
pairofperiodicorbitsforthedynamicalsystem,i.e.synchro- the oscillators. Thus, a non-zero value can be inferred even
nization is verified if ψ exists such that ψ = M(ψ ) and whenthereisnointeraction. Suchdiscrepanciescanbeover-
e e e
dM(ψ) <1. Theexistenceofthefixedpointψ isestab- come by careful surrogate testing [39, 40] – by rejection of
(cid:12) dψ |ψe(cid:12) e valuesbelowansurrogateacceptancethreshold,whichcanbe
l(cid:12)ishedthro(cid:12)ughthesimplealgorithmicproceduredescribedin
(cid:12) (cid:12) specifiede.g.asthemeanplustwostandarddeviationsamong
AppendixA.
manyrealizationofthemeasure.
In addition to the coupling strength and the directionality,
onecanalsoinferthecouplingfunctionthatcharacterizesthe
IV. DESCRIPTIONOFTHEINTERACTIONS interactions,i.e.thelawthatdescribesthefunctionalrelation-
ships between the oscillators. Its characteristic form reflects
Inferringthe parametersof the system notonlyallows for thenatureoftheoscillatorsandhowtheirdynamicsreactsto
evaluationofthesynchronizationasanepiphenomenoninits perturbations.
own right, but their probability X(c,Σ) also describes the The coupling function should be 2π-periodic. In the in-
N
interactionpropertiesoftheoscillators. Becausetheirdynam- ferentialframeworkunderstudy,thecouplingfunctionswere
ics is reconstructed separately, as described by Eq. (1), use decomposed into a finite number of Fourier components.
canbemadeonlyofthoseinferredparametersfromthebase The function describing the interactions between the two
6
(a) (b) (c) (d)
0
Figure 3: (Color online) Synchronization discrimination for the coupled phase oscillators (11). (a) Schematic of an Arnold tongue in the
coupling-frequencyε-ωplane: synchronizationexistsonlywithintheshadedarea[37]. (b)MapofM(ψ)forǫ12 =0.25demonstratingthat
theoscillatorsarenotsynchronized. (c)MapofM(ψ)foracasewherearootofM(ψ)=ψexists,i.e.wherethatthestateissynchronized.
(d)Thecorrespondingphasedifference,exhibitingtwophaseslips.
oscillators was decomposed by use of the odd parameters thetic signal, and parameter inference based on that signal)
q (φ ,φ ) c ,c ,... and the corresponding base func- the varianceof a particular parameterwould increase mono-
1 1 2 1 3
∈ { }
tionsΦ [q (φ ,φ )] sin(φ ,φ ),cos(φ ,φ ) uptoorder tonicallywithnoiseamplitude,asshowninFig.4. Thereare,
n 1 1 2 1 2 1 2
∈{ }
n of the decomposition. The reverse function q (φ ,φ ) however, a few notable exceptions. The inferentialcapabili-
2 1 2
∈
c ,c ,... wassimilarlydecomposed. ties rely on the volume of phase-space spanned by the vari-
2 4
{ }
ables. Astate ofsynchronizationwouldrepresentalimitcy-
cle fortheglobalsystem, andparameterinferenceof neither
oscillator would reach satisfactory precision. In such cases
V. TECHNICALASPECTSANDCONSIDERATIONS
a minimal amount of noise is typically needed, sufficient to
drivethesystemoutofequilibriumatleastonce. Duringthe
Thetechniqueisquitegenerallyapplicabletoabroadclass resultantphase-slip,thedatawouldbefillingthephasespace
of problems[54], and so there are a numberof technicalas- sufficientlyforcorrectparameterreconstruction.
pects and choices to bear in mind. We now discuss three of c. Time resolution. We nowsummarizethelimitsofan
them in particular: the number of base functions to be em- idealized data acquisition. The time step h is much bigger
ployedin theinferenceprocess(whichis partofthemodel);
theintensityofthenoisecharacterizingthesystem(whichis
anexternallyimposedconstraint);andthetimeresolution.
7 (a)
a. Number of base functions. Selection of the optimal
setofbasefunctionstodescribetheproblemisfarfromtriv-
ial. In general one wishes to have the minimal set that de-
scribes sufficiently well the model to be tested. Where the ω1
6.9
length of the data series is very long or effectively infinite,
onecanincludeanexcessivenumberofbasefunctionswith-
out immediate penalties. In reality, however, any unneeded
basefunctionjeopardizestheprecisionofthecoefficientsthat 0.001 0.1 0.2 0.3 0.5 1
E
really are relevant for the model, and the picture is further
1.8 (b)
complicatedwhenthemodeltobeadoptedis expectedto be
anoutcomeoftheinferencemachinery.Whereonedealswith
a longdata series, possiblywith a highsignal-to-noiseratio, 2 1.4
1
a relative large number of base functions can be used. The ε
speedofcomputationisalsoanimportantaspecttokeepun-
1
der consideration, given that having a large number of base
functionsvastlyincreasestheparameterspace,andthatitera-
1 0.001 0.1 0.2 0.3 0.5 1
tivecalculations(especiallymatrixinversion)slowthespeed
E
ofprocessingbythethirdpowerofthenumberofcoefficients.
Note that, even though the Bayesian inference is generally
Figure4: Statisticsoftheinferredfrequencyω1andcouplingε1for
popular in real-time applications, computational speed limi-
differentnoiseintensitiesE. Thesignaltobeanalysedisgenerated
tationsmean thatourinferenceframeworkforgeneralphase from Eq. (12). The dotted line shows the actual values of the pa-
dynamicscannotyetbeusedinthisway. rameters. The boxplots refer to the descriptive statistics (median,
b. Role of noise intensity. In general, the greater the quartiles, max. and min.) of 104 different runs of the generation-
noise intensity, the bigger the covariance of the inferred pa- and-inferenceloop.
rameters.Forarepeatedexperiment(e.g.generationofasyn-
7
7.1 (a) 0.5s VI. APPLICATIONS
5s
6.9 50s
ω(t)1 100s perTfhoermteacnhcneioqfuteheisafilgrostriathpmpl,ieadndtothseynnrtehaeltidcatdaaataretoantaelsytztehde.
6.7
Tocreatethesyntheticdata,weusedbothnumericalandana-
logueelectronicsimulations.
6.5

0 1000 2000
(b)
2
1.5
A. Numerically-generatedtestdata
ε(t)2 1
0.5 Numericallygenerateddatawereobtainedfrommodelsof
phaseoscillatorsandlimit-cycleoscillators.
0 1000 2000
Time [s]
Figure5:Inferenceof(a)atime-varyingfrequencyand(b)coupling
parameter from time series data generated by model (12) for four 1. Phaseoscillators
differentlengthsoftheinferencewindows.Thesizesofthewindows
areshowninthebox.
The phase oscillator model provides a sufficient basis for
the description of synchronization while being, at the same
time, analytically traceable. We thus test the detection of
synchronization (as explained in Sec. III) through Bayesian
than the correlationtime of the noise (assuming the latter to inferenceof synthetic data whose synchronizationis already
benon-zero). Atthesame time, however,h ismuchsmaller known. Themodelisgivenbytwocoupledphaseoscillators
thananyofthe sequentialtime-windowsusedasdata blocks subjecttowhitenoise
for inference, so that each block contains many data points.
Also, h ismuchsmaller than eitherof the oscillatorperiods. φ˙ =ω +ǫ sin(φ φ )+ξ (t) i,j =1,2. (11)
i i ji j i i
Each inference block is big enough to contain many cycles −
ofthedynamics(inparticular,morecyclesthatthosetypical
The parametersare ω = 1.2, ω = 0.8, ǫ = 0.1; param-
of a phase slip) while, at the same time, each block is small 1 2 21
eter ǫ is chosen so that the system lies close to the border
enoughtoprovidethedesiredresolutionofparameterchange. 12
of the Arnoldtongue (eitherjust inside or just outside). Be-
Itcanhappenthatthetimeresolutionofthechangeindy- cause we aim to demonstrate the precision of synchroniza-
namicalparametersisincompatiblewithanacquisitiontime- tion detection, we add no time-variability to the model, the
windowthatwouldguaranteeprecisionforotherparameters. inference is applied to a single block of data, and there is
The choice of the time-windowmust thereforebe done on a no spatial noise correlation with noise intensities E = 2.
ii
case-by-casebasis,dependingonthetypeofinformationthat The dynamics of the phase difference is described as ψ˙ =
is required from the system. Fig. 5 illustrates such a com- ∆ω ǫsin(ψ)+ξ (t)+ξ (t),where∆ω = ω ω isthe
1 2 2 1
promise. Weusethenumericalmodel(12)thatwillbeintro- frequ−encymismatch and ǫ = ǫ +ǫ is the res−ultant cou-
21 12
ducedinSec.VIA2toinvestigatethetime-resolutionforthe pling.Itisevidentthattheanalyticconditionforsynchroniza-
casewherethefrequencyω1(t) = ω1+A˜1sin(ω˜t)andcou- tion,i.e.theexistenceofastableequilibriumsolutionψ˙ <0,
plingamplitudeε (t) = ε +A˜ sin(ω˜t)werevaryingperi- is ∆ω/ǫ < 1. For ǫ = 0.25 (outside the Arnold tongue)
2 2 1 12
odicallyatthesametime.Theparameterswere:ω =2π1.1, the reconstructed map M(ψ) (Fig. 3(b)) after parameter in-
1
ω = 2π2.77, ε = 0, ε = 1, ω˜ = 2π0.002, A˜ = 0.1 ference has no root M(ψ ) = ψ : hence the oscillators are
2 1 2 1 e e
A˜ =0.5andnoisestrengthsE =E =0.15.Theparame- not synchronized. When ǫ = 0.35, even though the sys-
2 1 2 12
ters were reconstructed using four different window lengths temwasinsidetheArnoldtongue,noisetriggeredoccasional
for the inference. The results presented in Fig. 5 demon- phaseslips(see Fig.3(d)). We tested synchronizationdetec-
stratethat,forsmallwindows(0.5s),theparametersaresparse tiononthesame signalsusingthe methodsalreadyavailable
and sporadic, while for very large windows(100s) the time- intheliterature,basedonthestatisticsofthephasedifference
variability is faster than the size of the window and there is [5–7], but none of them was able to detect the presence of
cut-off on the form of the variability. The optimal window synchronizationundertheseconditions. Inspiteofthephase
lengthwillliebetweenthesetwo. Anotherinterestingfeature slips, our technique correctly detects the root M(ψ ) = ψ
e e
isthat,forthesmallestwindow(0.5s),thecouplingamplitude fromtheinferredparameters,revealingthattheoscillatorsare
improves with information propagation as time progresses, intrinsically synchronized as shown in Fig. 3(c): the phase
while the frequencyinferred(asa constantcomponentwith- slipsareattributablepurelytonoise(whoseinferredintensity
outbasefunction)remainssparsethroughoutthewholetime isgivenbythematrixE ),andnottodeterministicinterac-
i,j
interval. tionsbetweentheoscillators.
8
thereconstructedtime-variableparametersissatisfactorydur-
ingthenon-synchronizedintervals. Duringthesynchronized
intervals,however,theoscillatorsdonotspansufficientphase-
space to allow precise inferenceof the parameters(Fig. 6(a)
and (b), dashed lines). Within these synchronized intervals,
theposteriorprobabilitydistributionoftheparameterswasnot
peaked; however, it was sensibly differentfrom zero only in
that parameter region for which the correspondingnoiseless
dynamics is synchronized. Hence, despite the impossibility
of accurate parameter tracking, the detection of a synchro-
nized state (s(c) = 1) is always precise (Fig. 6(a) and (b),
greyshadedregions).
Figure6:(Coloronline)Extractionoftime-varyingparameters,syn- The reconstructed sine-like functions q (φ ,φ ) and
1 1 2
chronizationandcouplingfunctionsfromnumericaldatacreatedby q (φ ,φ ) are shown in Figs. 6(c) and (d) for the first and
2 1 2
(12). Thefrequency ω1(t)(a) and coupling ε2(t) (b) areindepen- secondoscillators,respectively. Theydescribethefunctional
dentlyvaried.Thedottedandfulllinesplottheparameterswhenthe
formoftheinteractionsbetweenthetwoPoincarésystemsin
twooscillatorsaresynchronizedforpartofthetime(ε1 =0.3),and Eq. (12). The reconstructed form of the coupling functions
notsynchronizedatall(ε1 =0.1),respectively. Theregionsofsyn-
wasevaluateddynamicallyforeachblock.
chronization,foundbycalculationofthesynchronizationindex,are
indicatedbythegrayshadedregions. (c)and(d)showthecoupling Next, the method was applied to deduce the predominant
functions q1(φ1,φ2) and q2(φ1,φ2) for time windows centered at direction of coupling as specified from the norm of the in-
t = 350s. Inbothcases,thewindowlengthwastw = 50sandthe ferred coupling base parameters. To illustrate the precision
couplingwasε12 =0.1. of the directionality detection, the frequencies were now set
constant, while both of the coupling strengths remained dis-
cretely time-varying. The parameters were ω = 2π1.3,
1
ω = 2π1.7, E = E = 0.2, andthecouplingfunctions
2 11 22
2. Limit-cycleoscillators
were,asinthepreviousexample,q (x ,x ,t) = x x and
i i j i j
−
q (y ,y ,t)=y y .Synchronizationwasnotreached,how-
i i j i j
To demonstrate the capabilities of the technique in trac- ever,forthesepa−rameters. Thecouplingsundergochangesat
ingtime-varyingparameters,couplingfunctions,directional- particular times, but otherwise remain constant, as shown in
ity and synchronization, we analyzed data from a numerical Fig. 7. The detected directionality index D was consistent
modeloftwocoupled,non-autonomous,Poincaréoscillators withtheactualvalues. Notethat,forunidirectionalcoupling,
subjecttowhitenoise, Ddoesnotquitereachunityonaccountofthenoise.
Tofurtherinvestigatetheabilitytotracksubtlechangesof
x˙ = r x ω (t)y +ε (t)q (x ,x ,t)+ξ (t),
i − i i− i i i i i j i time-varyingcouplingfunctions,weusedthesamemodelas
y˙i = riyi+ωi(t)xi+εi(t)qi(yi,yj,t)+ξi(t), (12) inEq. (12) to generatea syntheticsignalwherethe coupling
−
functionsareabsolutevaluesofthestatedifferencetoapower
r =( x2+y2 1) i,j =1,2.
i q i i − ofthetime-varyingparameter:
We tested several possibilities for the parameters: while let-
q (x ,x ,t)= (x x )ν(t) ,
ting the frequencies ωi and coupling parameters εi be time- x,i i j | j − i | (13)
varying, we ran numerical experiments with the coupling q (y ,y ,t)= (y y )ν(t) ,
y,i i j j i
| − |
functionq eitherfixedortime-varying.
i
As a first numerical experiment, we considered bidirec- where i,j = 1,2 and i = j. The exponent parameter
{ } 6
tional coupling (1 2), where the natural frequency of the
↔
first oscillator, and its coupling strength to the second one, 1
varyperiodicallyatthesametime:ω1(t)=ω1+A˜1sin(ω˜1t) ε =0.5
εand=ε20(.t1),ω= ε=2 +2πA˜12,sωin(=ω˜22t)π.1T.1h4e,oA˜the=r p0a.r2a,mA˜eter=s w0.e1r3e:, 0.5 ε21=0 εε1==00..24 εε21==00..33
2 1 2 1 2 D 0 2
ω˜ =2π0.002,ω˜ =2π0.0014andnoiseE =E =0.1.
1 2 11 22 ε =0.1
The coupling function was proportional to the difference in 1
−0.5 ε =0.5
thestatevariables: q (x ,x ,t) =x x andq (y ,y ,t) = 2
i i j i j i i j
−
yi−yj (thesamecouplingfunctionwasusedforconstruction −1 500 1000 1500 2000 2500
ofFig.1,Fig.4andFig.5). Thephaseswereestimatedasthe Time [s]
anglevariableφ = arctan(y /x ). Withε = 0.1,inastate
i i i 1
ofnosynchronization,thetime-varyingparametersω1(t)and Figure 7: Directionality of coupling D for discretely time-varying
ε2(t) are accurately traced as it can be seen in Fig. 6(a) and couplingamplitudesε1 andε2. Differentunidirectionallyandbidi-
(b). The coupling amplitude of ε1 = 0.3 corresponds to a rectionallycoupledstatesarereachedfordifferentvaluesofε1 and
state of intermittent synchronization, where the two oscilla- ε2,asindicatedbythesquareinsets.
tors are synchronized for part of the time. The precision of
9
Figure8: (Color online) Time-evolution of thecoupling function frommodel (12) withexponential timevariations (13). (a)-(d) Coupling
functionq1(φ1,φ2)fromthefirstoscillatorforfourconsecutivetimewindows. Thewindowlengthwastw =50s. Forsimplicityandclarity
onlythefunctionq1(φ1,φ2)isshown.Thebehaviorofq2(φ1,φ2)fromthesecondoscillatorwassimilar.
varied linearly with time ν(t) = 1 3 , and the other wereε = 0.7,ω = 2π15.9,ω = 2π17.5,A˜ = 0.03,ω˜ =
1 2 1
{ → }
parameters remained constant: ω = 2π1, ω = 2π2.14, 2π0.2andc=100isconstantresultingfromtheanaloguein-
1 2
ε = 0.2, ε = 0.3 and E = E = 0.05. The recon- tegration.Thephaseswereestimatedasφ =arctan(x˙ /x ).
1 2 11 22 i i i
structedphase couplingfunctionsqi(φ1,φ2) were calculated Forthegivenparameterstheoscillatorsweresynchronized,
from the inferred parametersfor the interacting terms of the sothattheseconddrivenoscillatorchangeditsfrequencyfrom
basefunctions. Theresultsforfourconsecutivewindowsare beingconstanttobeingtime-varying.Applyingtheinferential
presentedinFig.8. Itcanreadilybeseenthattheircomplex techniqueshowed, correctly, that the oscillatorswere indeed
form now is not constant, but varies with time. Comparing synchronized (s(c) = 1) throughout the whole time period.
theminneighboring(consecutive)pairs: (a)and(b),then(b) The frequencyof the driven oscillator was inferred as being
and(c),then(c)and(d),onecanfollowthetime-evolutionof time-varyingFig. 9(c). Performinga simple FFT (Fig. 9(d))
the functionalform. Even though we can follow their time- showed that ω (t) is periodic with period 0.2Hz (exactly as
2
variability,thetwomostdistantfunctionsFig.8(a)and(d)are set on the signal generator). Clearly, the technique reveals
of substantially differentshapes. Note also that, beside their information about the nature and the dynamics of the time-
form, the functions’ norm i.e. coupling strength also varies variability of the parameters – and is able to do so using a
(cf.theheightofthemaximainFig.8(a)and(d)). morerealisticsignalthanthatfromanumericalsimulation.
B. Analoguesimulations
C. Cardiorespiratoryinteractions
We also tested the technique on signals emanating from
Havingtestedourtechniqueontwoquitedifferentkindsof
analogmodels.Thesearereal,highlycontrollable,oscillatory
synthetic data, we now apply it to a realphysiologicalprob-
systemsandthenoiseontheirsignalsisrealratherthancon-
lem,toinvestigatethecardiorespiratoryinteraction.Theanal-
trived, as in the case of numerical models. It is attributable
ysis of physiological signals of this kind has already been
to environmental disturbances, thermal fluctuations, and the
foundusefulinrelationtoseveraldifferentdiseasesandphys-
inherentnonidealitiesof the circuit components. During the
iologicalstates(seee.g.[42]andreferencestherein). Transi-
process of data acquisition and discretization, measurement
noise can be introduced as well – noise which has no links
with the actualdynamicsof the interactingoscillators. Such
signalsprovideagoodtestofouranalysiscapabilities.
We analyzeddatafroman analogexperimentalsimulation
of two coupled van der Pol oscillators. Details of the elec-
tronic implementation are given elsewhere [41]. The noise
here arises mainly from the imperfections of the electronic
componentsandthereisalsomeasurementnoise.
Fig.9(a)showsthephaseportraitderivedfromthefirstos-
cillator,withtime-varyingfrequency,whichdrivesthesecond
oscillator
1x¨ µ (1 x2)1x˙ +[ω +ω˜ (t)]2x =0, Figure 9: Analysis of signals from an analogue simulation of the
c2 1− 1 − 1 c 1 1 1 1 system(14). (a)Phaseportraitfromtheoscilloscope;(b)frequency
c12x¨2−µ2(1−x22)1cx˙2+ω22x2+ε(x1−x2)=0, (14) ω˜1(t) from the external signal generator; (c) detected frequency
ω2(t) of the second driven oscillator; (d) Fast Fourier Transform
wheretheperiodictime-variabilityω˜ (t) = A˜ sin(ω˜t)(Fig. (FFT)ofthedetectedfrequencyω2(t).
1 1
9(b))comesfromanexternalsignalgenerator.Theparameters
10
2π
(a)
N
Ψ2:
0

1 (b)
I sync 2:9 2:8
0 0.22
(t)h0.9 0.2f(t)r
f 0.8 fh fr (c) 0.18

0 400 800 1200
Time [s] Figure11:(Coloronline)Couplingfunctionsinthecardiorespiratory
interaction calculated at different times. (a)-(c) Coupling function
Figure10: Synchronizationandtime-varyingparametersinthecar- q1(φ1,φ2)fromthefirstoscillator, and(d)-(f)q2(φ1,φ2)fromthe
diorespiratoryinteraction. (a)Standard2:N synchrogram. (b)Syn- secondoscillator.Thewindowtimeintervalswerecalculatedat:t=
chronization index for ratios 2:8 and 2:9 as indicated. (c) Time- 725sfor(a)and(d);t = 1200sfor(b)and(e);andatt = 1250s
evolutionofthecardiacfh(t)andrespiratoryfr(t)frequency. for(c)and(f).Thewindowlengthwastw =50s.
correspondingsynchrogramshowninFig.10(a),butprovided
tions in cardiorespiratorysynchronization have been studied
aclearerandlessambiguousindicationofsynchronization.
inrelationtoanæsthesia[43]andsleepcycles[44]. Itisalso
Thefunctionalrelationshipsthatdescribethecardiorespira-
knownthatmodulationsandtime-varyingsourcesarepresent,
toryinteractionsareshowninFig.11. Evaluatedforthreedif-
andthatthesecanaffectthesynchronizationbetweenbiolog-
ferenttime windows, the upperfigures(a)-(c)showthe cou-
ical oscillators [42, 45, 46]. For comprehensiveand reliable
analysisatechniqueisneededthatisable,notonlytoidentify plingfunctionq1(φ1,φ2)fromthecardiacoscillatingactivity,
thetime-varyinginformation,butwhichwillallowevaluation andthelower figures(d)-(f)show q2(φ1,φ2) fromthe respi-
rationoscillator. Theformofthefunctionsiscomplex,andit
of the interacting measures (e.g. synchronization and direc-
changesqualitatively over time – cf. Fig. 11(a) with (b) and
tionality), based solely on the information inferred from the
(c),or(d)with(e)and(f). Theinfluencefromrespirationto
signals. Wewillshowthatourtechniquemeetsthesecriteria.
heart (q (φ ,φ )) has a larger norm (i.e. coupling strength)
We analyse cardiorespiratory measurements from human 1 1 2
thanintheoppositedirection,indicatingthatthepredominant
subjectunderanæsthesia. Theirbreathingratewasheldcon-
directionofcouplingisfromrespirationtoheart.Onecanalso
stant,beingdeterminedbyarespirator. Forsuchsystemsthe
observethatq (φ ,φ )in(b)and(c)isofafairlyregularsi-
analytic model is unknown, in contrast to analogue and nu- 1 1 2
nusoidalformwith a stronginfluencefromrespiration. This
merical examples, but the oscillatory nature of the signal is
arisesfromthecontributionofthosebasefunctionsdescribing
immediately evident. The instantaneous cardiac phase was
thedirectinfluenceofrespiration(foradetaileddiscussionsee
estimatedbysynchrosqueezedwaveletdecomposition[47]of
[48]). Furthermore,Fig.11showsthatthefunctionalrelation-
theECGsignal.Similarly,therespiratoryphasewasextracted
ships for the interactionsof an open (biological)system can
fromthe respirationsignal. Thefinalphasetime-serieswere
inthemselvesbetime-varyingprocesses.
reachedafterprotophase-phasetransformation[12].
Application of the inferential technique reconstructs the
phase parameters that governthe interacting dynamics. Fig.
10(c) shows the time-evolution of the cardiac and respira- VII. GENERALIZATIONTONETWORKSOF
OSCILLATORS
tionfrequencies. Itis evidentthatthe constantpacingof the
breathingis well-inferred,andthat the instantaneouscardiac
frequency, i.e. “heart rate variability”, increases with time. Our parameter inference procedure can be applied with
Theinferredparameters,andtheircorrelations,isusedtode- only minimal modification to any number N of interacting
tect the occurrence of cardiorespiratory synchronizationand oscillatorswithinageneralcoupled-networkstructure.
thecorrespondingsynchronizationratio.Thesynchronization The notation of Eq. (1) is readily generalized for the N
evaluation I = s(c) 0,1 , shown in Fig. 10(b) re- oscillators, and the inference procedure, Eq. (7), is then ap-
sync
∈ { }
vealstheoccurrenceoftransitionsbetweenthesynchronized plied to the correspondingN-dimensionalphase observable.
and non-synchronizedstates, and transitions between differ- The number of base function should of course be refined to
ent synchronization ratios: from 2:8 (i.e. 1:4) at the begin- include terms relevant for the N-oscillator inference. Sub-
ning to 2:9 in the later intervals. Because the evaluation of procedures like the time-varying propagation, and the noise
the synchronization state is based on all of the given details inference,applyexactlyasbefore. Wenotethatthecomputa-
aboutthephasedynamics,theproposedmethodnotonlyde- tionalpowerrequiredincreasesexponentiallywithN,which
tectstheoccurrenceoftransitions,butalsodescribestheirin- makesthe methodunsuitablefor the inferenceof large-scale
herentnature. The results for I were consistent with the networks. However,forrelativelysmallnetworks,astandard
sync