WO2018087763A1 - Methods and systems for diagnostics - Google Patents

METHODS AND SYSTEMS FOR DIAGNOSTICS

BACKGROUND ART

References considered to be relevant as background to the presently disclosed subject matter are listed below:

[01] Daniels, B. C. (2005). Synchronization of Globally Coupled Nonlinear

Oscillators: the Rich Behavior of the Kuramoto Model.

[02] Goelman, G., R. Dan, F. Ruzicja, O. Bezdicek, E. Ruzicha, J. Roth and R. Jech (Submitted ). Frequency-phase analysis for resting-state fMRI. Human Brain Mapping

[03] Granger, C. W. J. (1969). "Investigating causal relations by econometric models and cross-spectral methods." Econometrica 37(3): 424-438.

[04] Lancaster, J. L., M. G. Woldorff, L. M. Parsons, M. Liotti, C. S. Freitas, L. Rainey, P. V. Kochunov, D. Nickerson, S. A. Mikiten and P. T. Fox (2000). "Automated Talairach atlas labels for functional brain mapping." Hum Brain Mapp 10(3): 120-131.

[05] Wu, X., X. Yu, L. Yao and R. Li (2014). "Bayesian network analysis revealed the connectivity difference of the default mode network from the resting-state to task- state." Front Comput Neurosci 8: 118.

[06] Panzica, F., Varotto, G., Rotondi, F., Spreafico, R. & Franceschetti, S. Identification of the Epileptogenic Zone from Stereo-EEG Signals: A Connectivity- Graph Theory Approach. Frontiers in neurology 4, 175, doi: 10.3389/fneur.2013.00175 (2013).

[07] Wang, M. Y. et al. Identification of the epileptogenic zone of temporal lobe epilepsy from stereo-electroencephalography signals: A phase transfer entropy and graph theory approach. Neurolmage. Clinical 16, 184-195, doi: 10.1016/j .nicl.2017.07.022 (2017). [08] Mao, J. W. et al. Dynamic Network Connectivity Analysis to Identify Epileptogenic Zones Based on Stereo-Electroencephalography. Frontiers in computational neuroscience 10, 113, doi: 10.3389/fncom.2016.00113 (2016).

[09] Goelman, G. & Dan, R. Multiple-region directed functional connectivity based on phase delays. Human brain mapping 38, 1374-1386, doi: 10.1002/hbm.23460 (2017).

[10] Goelman, G. et al. Frequency-phase analysis of resting-state functional MRI. Scientific reports 7, 43743, doi: 10.1038/srep43743 (2017). Acknowledgement of the above references herein is not to be inferred as meaning that these are in any way relevant to the patentability of the presently disclosed subject matter.

BACKGROUND

Understanding the nonlinear interaction between multiple coupled time-series functions (or in wider prospective multiple sequential series of values) have been proven essential to reveal and understand organization, function and interrelations between regions of various types of examined entities, such as nervous systems of subjects, seismological entities, and/or economical markets. Interaction between time- series that is commonly presented as networks and connectivity have been widely uses to understand brain organization and function . For example MRI data is used to construct structural (anatomical links), functional and effective connectivity between brain regions. The characterization of network connectivity is usually addressed in terms of the distinction between functional and effective connectivity. Functional connectivity typically refers to statistical dependencies amongst measured time-series, while effective connectivity typically rests on an explicit model of how those dependencies were caused (e.g., so called dynamic causal modelling and structural equation modelling).

Within the functional connectivity, there is a key distinction between directed and undirected functional connectivity. Simple measures such as correlations and coherence return undirected measures of statistical coupling, while procedures that appeal to temporal precedence afford inferences about directed functional connectivity. These include Granger causality (Granger 1969), transfer entropy and the procedure introduced here that is based upon significant phase delays.

Presently, network analysis is increasingly advancing the field of neuroimaging. Conventional construction of nervous system (e.g. neural network) models are generally based on statistical identification of pairwise interactions between respective pairs of regions of the nervous system, which it to be modeled, with an assumption of linear relations between them.

GENERAL DESCRIPTION

The present invention provides a novel technique for modeling the interaction between coupled multiple sequential-series functions including obtaining their hierarchal organization in various types of entities, such as seismological , economical and/or neuroscience entities, such as nervous systems. The technique of the present invention is in the field of high order statistics and enables to determine functional connectivity (linear as well as none linear connectivities) among multiple regions of the explored entity/system, as well as determining the directed functional connectivity between the regions, hierarchical/temporal order of the regions, and also for identifying region(s) sourcing a disorder in the system/entity (such as a seizes). Utilizing data indicative of a plurality of n signals (more than two signal components, typically four or more) obtained/acquired from a respective plurality of distinct regions of the explored system/entity (e.g. an individual's nervous system), the technique of the present invention determines statistical relations of higher statistical order between the plurality of more than two (typically 4 or more) distinct regions of the explored entity (e.g. nervous system).

According to some implementation of this technique specifically designed wavelet and spectral coherence analyses are used to characterize statistical dependencies between the n signal components of the nervous system (e.g. statistical dependencies over time) in terms of wavelet decompositions and systematic differences in phase relationships.

The technique of the present invention transcend the normal (linear) characterization of dependencies to consider high order moments. This allows to construct measures of functional connectivity among more than two nodes/regions of the explored system (e.g. a person's nervous system or muscle system) and assign an ordinal or hierarchical ranking to the implicit functional connectivity architecture.

The present invention, in its one prominent aspect, provides a method for use in identifying phenomenon/disorders in a subject/system, by extracting (non-invasively) marker/indicator obtained by processing the data indicative of the signals measured from the subject/system (e.g. a subject's muscle/nervous systems or parts thereof, or earths tectonic system, economical systems and/or other network behaving systems). The method includes:

(a) Providing data indicative of a plurality of n signals that were respectively acquired from n>3 distinct regions of the system/entity which is to be monitored/analyzed; and

(b) Processing at least one subset of m signals of the provided n signals, which are associated respectively with m the n regions (whereby n≥m and m≥3), to determine a connectivity network model of the m regions, The processing to determine the connectivity network model of the m regions includes the following:

i. Determining low order functional connectivity (FC) network model of said m regions, with order 0 = 1. This is achieved by processing pairs of the m signals (which are respectively associated with pairs of the m regions) to determine respective pairwise coherences of the respective pairs of signals. Each pairwise coherence of a pair of signals presents (or is indicative of) an

FC sub-network of order O = 1 in the m regions. Thus the pairwise coherences of the pairs of signals present an FC network model of the order 0 = 1 between the m regions; and

ii. Utilizing the FC network model of the order O between the m regions to determine a higher order FC network model of the m regions. This includes determining one or more coherent relationships of one or more of the FC sub-networks of the order O with either another one or more FC sub-network of order O or lower, and/or with one or more of the m signals. According to the invention, these coherent relationships are indicative of a higher order FC sub-networks of order O→ 0+1. To this end, the coherent relationships presenting a higher order FC network model of the m regions of order O→ 0+1. The method thus provides for obtaining a high order connectivity network model of the m regions. The model may serve as an indicator to the type/existence of the phenomenon in the inspected system/subject.

According to some embodiments of the present invention, the method further comprises determining a FC network model of a further higher order, of the connectivity between the m regions. This includes:

iii. repeating operation (ii) for one or more times whereby in each time utilizing the an FC network model obtained in the preceding operation (ii) to determine a further higher order FC network model; wherein determining the further higher order FC network model includes: determine respective coherent relationships of one or more of the higher order FC sub-networks obtained in the preceding operation (ii) with at least another one of the FC sub-networks or at least another one of the m signals. The coherent relationships between the higher order FC sub-networks and the another one of the FC sub-networks (of the same or low order) as well as the coherent relationships between the higher order FC sub-networks and one or more of the m signals, are indicative of a further higher order FC network model of the m regions.

According to some embodiments the method is configured to provides a complete FC network model of the m regions. The complete FC network model being a model having a single sub-network of size L that equals the number of regions m. This In this regards, in case m = 3 or m=4 the higher order FC network model obtained in (ii) presents the complete FC network of the m regions. Otherwise, in (iii) the method includes repeating of the operation (ii) for 0-1 times, whereby O being an order of the complete FC network model of the m regions (O is an integer satisfying O > Log₂(m)).

It should be understood that according to some embodiments operations (i) and the one or more repetitions of (ii) are performed in combination (e.g. in a combined processing of the m signals).

It should also be understood that the method of the present invention may be implemented non-invasively by processing the signals acquired externally from the system.

According to some embodiments the method is conducted on time window portions of the signals. In specific implementations, it may be preferably that the connectivity network ids determined for simultaneous time window portions of the m signals (namely associated with the same timing). In such embodiments, the method includes extracting from the m signals, m respective time window portions associated with corresponding time slot, and utilizing the m respective time window portions to carry out the processing of the subset of m signals. In some implementations the time windows are selected based on time slot data indicative of at least one of the following: a timing of the time slot at which to extract the corresponding time window portions of the m signal, a duration of the time slot, and a resolution of the signal portions to be extracted. The time slot data may be obtained by at least one of the following:

(i) provision of event timing data indicative of a timing of occurrence of event associated with said phenomenon and use of the event timing data to set the timing of the time slot;

(ii) setting the duration and/or the resolution of the time slot based on a type of said phenomenon; and

(iii) receiving input data indicative of at least one of said: (a) the timing, (b) the duration, and/or (c) the resolution.

According to some embodiments determining of the functional connectivities between m regions includes applying time-frequency transformations to at least time window portions of the m signals to obtain m respective transformed representations being respectively indicative of frequency contents of at least said time window portions of the m signals in at least one frequency band. In specific examples the time- frequency transformation include wavelet transform, and the transformed representations are provided as complex wavelet representations. In some embodiments the at least one frequency band is a predetermined frequency band. In some embodiments the at least one frequency band is determined based on a type of the phenomenon/disorder that is to be characterized. . In some embodiments the at least one frequency band is provided as input.

According to some embodiments of the present invention the method includes identifying the phenomenon/disorder by comparing the connectivity network model of the m regions with reference data indicative of expression of functional networks between the m regions in similar systems having such phenomenon/disorder. Additionally the method may further include identifying the phenomenon/disorder by comparing the connectivity network model of the m regions with reference data indicative of expression of functional networks between the m regions in similar systems not having the phenomenon/disorder. It should be noted that the above discussed connectivity network model of the system, may be functional, hierarchical and/or directed connectivity model.

To this end some embodiments of the present invention are further exploit the above described invented technique of determining functional connectivity (FC) between regions, in order to determine the hierarchical order of the regions. This is based on the inventor's understanding that the different regions of a network behaving system, such the nervous system, behave as a network in the sense that they source coupled time-signals (e.g. the EEG signals), whose hierarchical temporal order (indicative of the of the hierarchical order of the regions), can be determined by the nonlinear coherence (i.e., phases) between frequency components of the signals sourced from different regions in one or more frequency bands (see for example Ref [09]).

Thus, according to some embodiments of the present invention the method further includes carrying the following operation in order to determine the hierarchical connective network model of the system:

(c) determining temporal hierarchy in the connectivity network between said m regions by computing relative phases between (m-l) regions of the m regions relative to another region of the m regions and wherein the relative phases are indicative of the temporal hierarchy between the m regions.

According to some embodiments of the present invention, the method further includes carrying the following operation in order to determine the directed connective network model of the system: providing assumption data indicative of temporal hierarchical relative order of two of the m regions and processing the assumption data together with said phase relations of the (m-l) regions to determine directional hierarchy of the connectivity network thereby modeling directed functional connectivity network of said m regions.

According to some embodiments of the present invention, the method further includes determining a type of said connectivity network indicative of whether said network is "linear", "non-linear" or a "disconnected" network. Here a "linear" network is defined as a 1 ^st order network composed of multiple pairwise coherences where each might be related to same or different time-frequency windows. The above described technique for identifying disorders in systems such as a nervous system (e.g. the brain), a muscle system (e.g. the heart system), or other network behaving systems, is further described in more detailes in the detailed description chapter below. The connectivity network presents a clinical indicator indicative of a phenomenon being a nervous system disorder.

The inventor's of the present invention have further developed this technique of the present invention for identification of a region sourcing the disorder (e.g. a region of source of the diorder such as a region sourcing a seizure). The method is based on the assumptions that:

(1) some disorders in network behaving systems, which may be characterized as seizure like disorder, start/sourced at specific location(s) and spread like waves theough the functionally connected regions. This gives rise to a so called global temporal 'first' region being the source of the diorder; and

(2) there are no global temporal 'last' (i.e., specific region(s) in which seizures/diorders stop).

To this end, as indicated above the present invention by processing a provides a technique for determining the functional connectivity and hierarchy (e.g. temporal) between a subset of m>2 signals (e.g. 3 siganls, 4 signals or more) measured from m different regions of a system. Presently the method is implemented efficiently (with reasonable computational complexity) for determining the functional connectivity as well as the hierarchy of subsets of m=3 regions/signals or m=4 regions/signals or also the functional connectivity and hierarchy of larger subsets. As described above, the functional connectivity in the subset of m signals/regions and their hierarchy may is determined based on nonlinear coherences between the multiple m signals. In this regards the inventors of the present inventiomn have developed a novel technique for harnesting the information about coherent relashionships between m signals of a signal subset, for determining the region of source (e.g. siezure source) of a disorder out of n>m regions from which signal are measured. This is achieved by processing the specifically selected plurality k of different subsets of m signals out of the total n measured, determining the coherent relashionships between the signals of each subset, and applying a certain statistical processing to the coherent relashionships obtained for the k subsets to determine one of the n regions that is the "global first", namely being the source of the disorder/siezuer. This provides for combining high order statistics to identify neural regions (i.e. electrodes) whose time signal precedes the others (i.e., identify the leading electrode/signal).

In some implementations, the technique of determining the source region of a disorder was implemented for processing EEG signal(s) obtained/measured from n neural regions of the brain (e.g. about n~100 regions) in order to identify a candidate region(s) whose signal precede most of other signals that are measured from other regions. In this case the candidate region(s) were considered to be the source zone/region of the neural disorder.

For example, as is further described and examplified below, the technique of the present invention was tested for processing EEG signals to determine/identify epileptogenic zones (regions in the brain of a subject) sourcing an epileptic seizure or its onset. Epilepsy affects 0.5-1% of the population worldwide. Estimates indicate that 20 - 30%) of patients with epilepsy are refractory to all forms of drug therapy. Some of these medically intractable patients are candidates for surgical treatment in an attempt to achieve seizure freedom. In the US alone 100,000 to 200,000 patients are estimated to be candidates for epilepsy surgery, with only about 4000 surgeries performed yearly. The goal of epilepsy surgery is to remove the abnormal area of cortex from which the seizures originate (i.e. epileptogenic zone) without causing any significant functional impairment. To identify this area (region of the brain), information from scalp EEG electrodes and from subdural and intracerebral depth electrodes (invasive EEG - iEEG) is typically used, in conjunction with anatomical (MRI) and functional (SPECT and PET) imaging. However, in many cases surgery fails to achieve seizure freedom. This is mainly because conventional techniques often fail to correctly identify the epileptogenic area (region of the brain) sourcing the epileptic seizure. Traditionally, the electrophysiological recordings (e.g. from EEG, iEEG and/or functional imaging) are interpreted by visual inspection of experienced epileptologists or physiologists. Alternaively in the recent years several algorithms were developed to identify spike- rate, seizures and seizure onset using a variety of approaches such as pattern recognition or graph analysis (see for example references [06] to [08] above). However, these conventional techniques often fail to correctly identify the region of the brain sourcing the epileptic disorder. This in turn leads to failure of many epileptic surgeries. Therefore, there is a critical clinical need for a reliable technique to aid in the identification of epileptogenic zone sourcing the epileptic disorder. Thus, according to some aspects/embodiments of the present invention, there is provided a novel technique for determining a region of sourcing a disorder in the system. The system may be any netwrok behaving system, of which different regions can be considered to be connected in a network, for example: the nervous system (e.g. the brain), and/or the muscle system (e.g. the heart) or other systems (e.g. tactonic system and/or economical system). For that matter the disorder whose source is to be identified may be for example an epileptic disorder, heart disorder, and/or other disorder such as earthquakes.

For example, according to some embodiments the method of the present invention is adapted for identifying a certain region of said n regions being a source region of the phenomenon. In such embodiments the method further includes the following operations:

(d) selecting from the n signals a group of m-1 signals whereby one of said m-1 signals of the group is selected as a reference signal;

(e) carrying out operations (b) and (c) described above n times, with n different respective subsets of m signals each including the group of m-1 signals selected in (d) and a different one of said n signals set as a candidate signal for being a source of said phenomenon. This provides for determining n relative phases between each one of the n signals being set as the candidate signal and the reference signal of the group of m-1 signals respectively;

(f) repeating operation (d) and (e) a plurality of k times with k different groups of m-1 signal, to determine per each candidate signal of said n signals, k relative phase between the candidate signal and reference signals of the k different groups. And further summing the k relative phases of the candidate signal to obtain an average relative phase of the candidate signal. This provides for obtaining n average relative phases each associated with a different one candidate signal of the n signals; and

(g) identifying the candidate signal(s) having the highest average relative phase as being associated with the source regions of the phenomenon/disorder in the system. In some implementations the technique for determining the source of the disorder is implemented by carrying out the above operations for several frequency bands and/or several time windows. To this end the method may further include: (h) carrying out operations (b) to (f) for each of a plurality J of selected frequency bands and wherein said summing of the k relative phases includes summing together absolute values of the J*k relative phases obtained between the candidate signal and the reference signals of said k different groups, for each of the J frequency bands to thereby obtain said average relative phase of the candidate signal computed for the plurality of frequency bands; and wherein said (g) identifying the candidate signal includes selecting one or more of the candidate signals having the highest average relative phase for the plurality of frequency bands as most likely candidates for being a source of said phenomenon. According to some embodiments in (e), the operations (b) and (c) are carried out n times at each time with a different candidate signal of said n signals, and with the similar group of m-1 signals as used for other candidate signals. The similar group of m-1 signals is characterized by the identity of the (m-1) signals and by the identity of the reference signal selected among the m-1 signals. Carrying out the method this way may provide that in the summation in (f) of the k relative phases of each of the candidate signal of the n signals, the average relative phases of candidate signals which are not associated with the source region of the phenomenon are suppressed.

In some embodiments the phenomenon is an epileptic disorder and the source of the phenomenon is a region sourcing an epileptic seizure in the system. According to some implementations of the invention, the technique described above for determining a the source region sourcing a disorder is used for the interpretation of intra- and extracranial electroencephalography (EEG). The technique can be used online (e.g. in real time) for intraoperative use e.g. during surgery and/or offline (e.g. in data post processing mode), to identify epileptogenic zones for epileptic patients who are candidates for surgical treatment.

To this end, differently from conventional EEG signal analyses tecnques, which are generally aimed at identifying changes in activity patterns indicated by the EEG signals, the technique of the present invention provides for determining a temporal hierarchy among the measured time-signals. As indicated above, the temporal hierarchy among the measured EEG time-signals is determined based on nonlinear coherences between multiple electrophysiological recordings. This provides for combining high order statistics to identify neural regions (i.e. electrodes) whose time signal precedes the others (i.e., identify the leading electrode). The neural regions, from which EEG signal(s) that precede most of other signals are measured, are considered to be the source zone/region of the neural disorder (e.g. presenting the epileptogenic zones sourcing the seizure or its onset).

According to another prominent aspect of the present invention there is provided a method for identifying a certain region of n regions of a system as a source region of a phenomenon in the system. The method includes:

(a) selecting from n signals obtained from n regions of the system a group of m-1 signals whereby one of said m-1 signals of the group is selected as a reference signal;

(b) determining high order connectivity networks between with n different respective subsets of the m signals, whereby each subset includes the group of m-1 signals selected in (a) and a different one of the n signals set as a candidate signal for being a source of the phenomenon and thereby determining n relative phases between each one of the n signals being set as the candidate signal and the reference signal of the group of m-1 signals respectively;

(c) repeating operation (a) and (b) a plurality of k times with k different groups of m-1 signal, to determine per each candidate signal of said n signals, k relative phase between the candidate signal and reference signals of said k different groups; and summing the k relative phases of the candidate signal to obtain an average relative phase of the candidate signal; thereby obtaining n average relative phases each associated with a different one of said n signals serving as candidate signal;

(d) identifying the candidate signal having the highest average relative phase as being associated with the source regions of said phenomenon.

According to yet another prominent aspect of the present invention there is provided a method for identifying a phenomenon in a system. The method includes:

(a) Providing data indicative of a plurality of n signals being respectively acquired from n>4 distinct regions of the system;

(b) Processing at least one subset of m signals of said n signals associated with m regions of said n regions, whereby n≥m and m≥4, to determine a connectivity network model of the m regions , whereby said processing to determine the connectivity network model of the m regions comprises: i. Determining low order functional connectivity (FC) network model of said m regions, with order O = 1, by processing pairs of said m signals, which are associated respectively with pairs of the m regions, to determine respective pairwise coherences of the respective pairs of signals, whereby the pairwise coherences of the pairs of signals are indicative of the FC sub-networks of order 0 = 1 between the m regions presenting an FC network model of the order 0 = 1 between the m regions; and

ii. utilizing said FC network model of the order O between the m regions to determine a higher order FC network model of order O→ 0+1, whereby said utilizing comprises processing pairs of the FC subnetworks of order O and determining respective coherent relationships between the FC sub-networks of each pairs of FC sub-networks , whereby the coherent relationships between the respective FC subnetworks of the pairs of FC sub-networks of order O are indicative of a higher order FC sub-networks of order O→ 0+1, presenting a higher order FC network model of said m regions;

thereby obtaining a high order connectivity network model of said m regions usable as an indicator of the phenomenon.

According to yet another prominent aspect of the present invention there is provided a system for identifying phenomenon. The system includes:

(1) a receiver module adapted for receiving data indicative of n signals acquired from n>4 distinct regions of an entity to be inspected; and

(2) a functional connectivity mapping module configured for processing at least one subset of m signals of said n signals associated with m regions of said n regions, whereby n≥m and m≥4, to determine a connectivity network model of the m regions. The functional connectivity mapping module is configured and operable to carrying out the following:

i. Determine low order functional connectivity (FC)network model of the m regions, with order 0 = 1. This is achieved by processing pairs of the m signals, which are associated respectively with pairs of the m regions, to determine respective pairwise coherences of the respective pairs of signals. The pairwise coherences of the pairs of signals are indicative of the FC sub-networks of order 0 = 1 between the m regions. These in turn present an FC network model of the order 0 = 1 between the m regions; and

ii. Utilizing said FC network model of the order O between the m regions to determine a higher order FC network model of the m regions. , This includes determining one or more coherent relationships of one or more of the FC sub -networks of order O with another one or more of the FC sub-networks (e.g. of order O or lower) or with another one or more of the m signals. Such coherent relationships are indicative of a higher order FC sub-networks of order O→ 0+1, presenting a higher order

FC network model of said m regions;

The system thereby provides a high order connectivity network model of the m regions.

According to some embodiments, the system further includes an identifier module capable of identifying said phenomenon by comparing said higher order connectivity network model with reference data indicative of expression of functional networks between said m regions in entities expressing said phenomenon.

According to some embodiments, the system further includes a directed connectivity mapping module capable of determining temporal hierarchy in said connectivity network model of said m regions. In this case the connectivity mapping module is configured and operable processing said higher order FC network model to determine phase relations of (m-1) regions of the m regions relative to another region of the m regions and utilizing said phase relations for determining the temporal hierarchy between the m regions.

In some embodiments the system further includes a classifier module capable of determining a type of the connectivity network indicative of whether the network is linear, non-linear or a disconnected network.

In some embodiments the system includes a Frequency-Time domain signal processor Module adapted for applying time-frequency transformations to at least time window portions of the m signals to obtain m respective transformed representations of the m signals. The transformed representations are indicative of respectively frequency contents of at least said time window portions of the m signals in at least one frequency band. The processing of the m signals by the functional connectivity mapping module is carried out on said m respective transformed representations of the m signals. In some implementations the Frequency-Time domain signal processor Module includes a wavelet transform module adapted to apply wavelet transform to said m signals to obtain said transformed representations as complex wavelet representations. In some implementations the , wherein said Frequency-Time domain signal processor Module the a TimeSlot-FrequencyBand provider module adapted to obtain data indicative of at least one time slot and at least one frequency band which are of interest for of interest for identifying said phenomena. In such cases the time-frequency transformations are conducted to extract portions of the m-signals which correspond to said at least one time slot and at least one frequency band.

According to some embodiments of the present invention the system includes a Disorder Source Identification module. The Disorder Source Identification module includes:

A. a subset generator module adapted for:

- selecting from said n signals, plurality of k different groups of m-1 signals and selecting one signal of the m-1 signals of each group as a reference signal;

- constructing a plurality of k*n different subsets of m signals each, whereby each subset includes one of said n signals selected as candidate signal and one of said groups of m-1 signals; and

- operating said functional connectivity mapping module and said directed connectivity mapping modules to determine hierarchical connectivity models per each one of said of k*n different subsets, whereby the hierarchical connectivity models are indicative of the respective K*n relative phases between the candidate signal and the reference signal of each of the n*K subsets;

B. a relative phase averager module configured and operable for summing the K relative phases associated with the same candidate signal of said n signals to obtain n averaged relative phases associated with said n signals respectively; and C. a candidate source region selector module configured and operable for selecting one or more of said n signals having the highest averaged relative phases as most likely candidates for being a source of said phenomena.

According to yet further another prominent aspect of the present invention there is provided a system for identifying a source region of a phenomena in an entity. The system includes:

(1) a receiver module adapted for receiving data indicative of n signals acquired from n>4 distinct regions of an entity to be inspected; and

(2) a Disorder Source Identification module comprising:

A. a subset generator module adapted for:

- selecting from said n signals, plurality of k different groups of m-1 signals and selecting one signal of the m-1 signals of each group as a reference signal;

- constructing a plurality of k*n different subsets of m signals each, whereby each subset includes one of said n signals selected as candidate signal and one of said groups of m-1 signals; and

- process the k*n different subsets to determine data indicative of the respective K*n relative phases between the candidate signal and the reference signal of each of the n*K subsets;

B. a relative phase averager module configured and operable for summing the

K relative phases associated with the same candidate signal of said n signals to obtain n averaged relative phases associated with said n signals respectively; and

C. a candidate source region selector module configured and operable for selecting one or more of said n signals having the highest averaged relative phases as most likely candidates for being a source of said phenomena.

According to further yet another prominent aspect of the present invention there is provided a system for identifying phenomenon in a certain entity (system/subject). The system includes:

(1) a receiver module adapted for receiving data indicative of n signals acquired from n>4 distinct regions of an entity to be inspected; and

(2) a functional connectivity mapping module configured for processing at least one subset of m signals of said n signals associated with m regions of the n regions, whereby ri≥m and m≥4, to determine a connectivity network model of the m regions by carrying out the following:

i. Processing pairs of said m signals, which are associated respectively with pairs of the m regions, to determine low order functional connectivity (FC) network model of order 0 = 1 between the m regions.

The processing includes determining pairwise coherences of the pairs of signals, being respectively indicative of FC sub-networks of size L = 2⁰⁼¹ between the m regions. The sub-networks presenting together an FC network model of the order 0 = 1 between the m regions; and ii. Processing the FC network model of the order O between the m regions to determine a higher order FC network model of order O→ 0+1. The processing includes determining pairwise coherent relationships of pairs of the FC sub-networks of the order O. The pairwise coherent relationships of the pairs of FC sub-networks of order O are indicative of higher order FC sub-networks of order O→ 0+1. The higher order FC sub-networks thereby presenting a higher order FC network model of order O→ 0+1 between said m regions;

thereby providing a high order connectivity network model of the m regions.

BRIEF DESCRIPTION OF THE DRAWINGS

In order to better understand the subject matter that is disclosed herein and to exemplify how it may be carried out in practice, embodiments will now be described, by way of non-limiting example only, with reference to the accompanying drawings, in which:

Fig. 1A is a flow chart schematically illustrating a method 100A for processing time series signals obtained from a plurality of regions of a system/entity to determine a functional connectivity network between the regions;

Fig. IB is a flow chart schematically illustrating a method 100B for processing a functional connectivity network to determine additional properties of the network such as hierarchical/directed connectivity thereof and linearity thereof, and/or identifying phenomenon indicators; Fig. 1C is a block diagram of a system for identification of an expression of a disorder/phenomena in a system/subject based on a model of connectivity network of the subject;

Figs. ID to IF are schematic illustrations of a functional connectivity networks obtained according to the technique of the present invention;

Figs. 1G and 1H are schematic illustrations of respectively hierarchal and effective connectivity networks obtained according to the technique of the present invention;

Figs. II to IK are schematic illustrations of several reference models of connectivity networks indicative of expression of normal and abnormal condition/phenomena in subjects;

Fig. 2 shows wavelet time-frequency contour-plots of four cortical regions of a subject and the functional connectivity (FC) between them determined according to the technique of the present invention;

Figs. 3A and 3B are a flowchart 100C and a block diagram 200B of a method and a system according to an embodiment to the present invention for identifying a region of source of a phenomena in a system;

Fig. 4 is a schematic illustration of the implementation of method 100C;

Figs. 5 and 6A to 6C are graphical illustrations of test study conducted using the technique of the present invention for identifying epileptogenic zone of a patient;

Fig. 7 is a graphical illustration_showing results of computer simulations of the technique of the present invention carried out by using the Kuramoto model of a coupled oscillators system to derive coupled time-series functions as a function of the coupling strength and noise;

Fig. 8 is a graphical illustration_showing dependency on Gaussian noise of the 4- region-FC and pairwise coherence phases;

Fig. 9 shows amplitude 3-seed-SPM versus joint-seed-SPM oi the ventral visual system for wavelet scale 1 (0.02-0.03Hz);

Fig. 10 shows a phase 3-seed-SPM oi the ventral visual system for wavelet scale 1 (0.02-0.03Hz).

Fig. 11 shows amplitude 3-seed-SPM versus joint-seed-SPM of the motor system for wavelet scale 1 (0.02-0.03Hz). Fig. 12 shows a phase 3-seed-SPM of the motor system for wavelet scale 1 (0.02-0.03Hz).

DETAILED DESCRIPTION OF EMBODIMENTS

The present invention provides a novel technique for modeling data in the form of sequential series functions (such as time-series) obtained from monitoring of a plurality of regions/measures (such as blood samples) of a certain system (e.g. an being a monitored object or a subject/"body"/entity) to identify/reveal/model nonlinear coupling/connectivity/ networks indicative of functional link between the regions/measures. The model of the coupling between the monitored regions of the monitored entity can be further used to identify/determine/diagnose various phenomena occurring/existing/characterizing the monitored entity/body.

In general, the technique of the present invention may be used in various disciplines, such as seismology, economy, medicine, data networks and/or neuroscience, in which modeling of coupled sequential/ time-series data obtained from multiple regions/measures may be used to identify certain phenomena/disorder/condition in the particular discipline. In this regards, the phrases region and subject! entity used herein should be interpreted broadly, to fit the discipline for which the technique of the present invention is employed. For example, in connection with the none limiting list of disciplines indicated above the phrase regions may designate as follows: geographical/seismologically-defined regions in seismological implementations of the present invention; markets/indexes/stocks in economical implementations of the present invention; blood test values in cancer patients; monitoring values in ER; and/or nervous system regions/nodes in neurology and neuroscience implementations. Accordingly the phrases body/entity may designate as follows: a planet or geographical zone thereof in seismological implementations; a markets in economical implementations; and/or a nervous system or part thereof in neurology and neuroscience implementations.

Accordingly, it should be understood that the phrases signal/data-series are used herein interchangeably to indicate data which is obtained/acquired from the monitored regions during a certain sequential order, and which indicates a certain ordered series of monitored property acquired from the monitored regions. Data-series may be construed herein as series of data-points taken in a sequential order (e.g. time) that are used in statistics, signal processing, pattern recognition, econometrics, mathematical finance, weather forecasting, earthquake prediction, neuroscience, imaging, astronomy, communications engineering, and largely in any domain of applied science and engineering which involves temporal and or sequential measurements.

There are various conventional techniques for analyzing a sequential series such as a time-series. These may be divided into frequency and time-domain methods. However, for many applications analysis of each time-series alone (i.e. of single time series) is not sufficient and knowledge on the interactions between time series may be essential. This is particularly important considering that in most cases time-series are coupled. The assumption of coupling is almost always true although the level of coupling varies between the different fields. Such knowledge has the potential for example to aid in understanding of changes in the stoke market, to predict earthquakes by studying coupled seismological waves, to predict cyclones and hurricanes by studying ocean waves, to define cardiac status by analyzing coupled cardiac waves or to identify onset location in epileptic seizures. Typically, conventional methods to analyze interaction between time-series are limited to interaction between pairs. However, even when all possible pairs, in an N coupled time-series system, are calculated, it is not sufficient to fully describe the interaction of the system. The subgroup of systems whose pair-interaction is sufficient is the 'linear' sub-group, but most systems of interest are nonlinear in this respect.

In the description below a coupled signals (e.g. coupled sequential/ time-series) obtained from a system/entity (e.g. from different respective regions thereof are considered as a network. Typically each series/signal presents, or is acquired from, a network node/region and network connectivity is given by the strength of coupling between the signals/series. Accordingly the phrase connectivity and functional connectivities used herein should be understood in the broad sense as coupling between sequential series. In this regards, it should be noted that in the following the term signals, and/or time series are used interchangeable to refer to sequential series, and accordingly these terms should be interpreted broadly to cover any type of sequential- series/signals and not only timely ordered series.

Furthermore, the directions of information flow correspond to the phases between series and therefore to the system sequential hierarchy. The characterization of network connectivity is usually addressed in terms of the distinction between functional and directed (effective) functional connectivity. Functional connectivity refers to statistical dependencies (coupling strength) amongst measured signals/time-series, while directed functional connectivity rests on causality; e.g. providing an explicit model of how those dependencies between the signals/time-series are caused (e.g., dynamic causal modelling, structural equation modelling or directed phase lag index). Within the functional connectivity, there is a key distinction between directed and undirected functional connectivity. Simple measures such as correlations and coherence return undirected measures of statistical coupling, while procedures that appeal to temporal precedence afford inferences about directed functional connectivity. These include Granger causality, transfer entropy and the procedure introduced here that is based upon significant phase delays.

In seismological implementations the sequential series data/signals may be acquired from sensors, such as the Global Seismographic Network and/or from other seismographs coupled for sensing the respective geographical regions. In economy, the sequential series data/signals may be for example exchange rates. In neuroscience, sequential series data/signals may be data/signals obtained from monitored regions of a nervous system by various techniques such as resting-state functional MRI (fMRI), stimulus-driven fMRI, electrophysiology, Magnetoencephalography (MEG) and/or Electroencephalogram (EEG).

The novel technique of the present invention provides for processing and/or transforming signals (sequential series data) measured/acquired from the plurality of monitored regions of the entity being monitored, to determine and obtain a connectivity network model (herein after referred to also just as model) indicative of functional and/or directed functional connectivity between the monitored regions.

In this connection it should be understood that the phrase functional connectivity refers to the mare statistical association between the activities of the monitored regions. The phrase directed functional connectivity refers to the directionality, i.e., the ordinal or hierarchical ranking, of the time-series between the monitored regions for which functional connectivity exists.

In various embodiments of the present invention the model of the statistical association (e.g. the model of the directed functional connectivities) between the monitored regions is further utilized in order to identify indicators/markers (e.g. patterns in the model) which may be indicative of one or more phenomena/disorders in the monitored entity/body and/or of the condition of the monitored subject. Accordingly, for example by comparing the connectivity network model obtained by the technique of the present invention with a reference model/data indicative of marker of the existence and/or inexistence of predetermined phenomena/disorder/condition, the phenomena/disorder/condition can be diagnosed, and/or determined to be existing or not, in the monitored subject/entity.

In this connection, advantageously, the present invention provides a novel method and system which are configured and operable for processing the time/ sequential series data/signals acquired from the plurality of regions to identify high order statistical association (e.g. statistical dependency/correlation/coherence) between the activities of these regions to orders higher than the 1^st statistical order. The high order (2^nd statistical order or higher) signal processing/transformation, which is used according to the present invention, will be explained in further details below. The use of such high order signal processing/transformation enables to convert the time- series/signals into a network connectivity model which reveals functional/hi erarchical/directed connectivities, that could, or would, not be revealed when relying on conventional techniques for determining functional/directed network connectivities based on only a first order statistical relations. Accordingly, the technique of the present invention enables to model the network connectivity of a subject/entity with improved details and/or with higher accuracy, and thereby facilitates more reliable and accurate identification and diagnostics of phenomena/disorder affecting the subject/entity and/or the condition of the subject.

Reference is made together to Figs. 1A to IK which illustrated and exemplify the technique of the present invention. Figs. 1A and IB are flowcharts a block diagram, schematically illustrating a methods 100A and 100B for determining existence and/or behavior of functional (100A) and directed (100B) connectivity networks between a plurality of regions of the monitored entity, and possibly also identifying a disorder/phenomena in the monitored entity. Fig. 1C is a block diagram exemplifying the configuration of a system 200A configured according to embodiments of the present invention carrying out the methods 100A and/or 100B, e.g. for use in identification of disorders/phenomena in a monitored subject/entity E. Figs. ID to IF are schematic illustrations of a modeled connectivity networks of a monitored subject/entity E. More specifically Figs. ID to IF illustrate the functional connectivities between plurality of monitored regions {Ri} of the subject/entity E, whereby Fig. ID illustrates the connectivity network obtainable via 1^st order (0=1) statistical processing/transformation of the time series data/signals obtained from the regions. Figs. IE and IF illustrate 2^nd and 3^rd order statistical processing/transformation carried out according to the present invention on the time series data/signals obtained from the regions, to determine/identify functional and/or directed connectivities, which might not have being revealed by the 1^st order statistical processing/transformation.

In this regards, it should be noted that the term order O used herein, is related to (indicative of) the maximal size of the connectivity network (or sub-networks) that can be characterized/identified by the processing. In the terminology used herein, connectivity network model of order O should therefore be understood as a model of m regions, which includes one or more connectivity networks (e.g. which are also referred to herein as subs-networks in case of a plurality) each being of order O or lower. As described in more details above and below, a connectivity network/subnetwork of order O may characterize the connectivity of up to L=2° or less. Thus to obtain a complete connectivity network of m regions (namely where the collective connectivity between all regions is characterized) statistical (coherence) processing of order 0=Log₂(m) should be carried out. In such cases the model will include a single connected network. To this end, according to the present invention, processing/computation carried out to order 0=1 is provides for characterizing pairwise connections between pairs of regions/signals. These each generally present connectivity networks of order 1 (linear network) while together they may be considered a connectivity network model order 1 of the m modeled regions/signals, that is being a linear model. However, according to the present invention, further iterations of higher order high order statistical computation is carried out for determining non-linear (higher order) connectivity networks/subnets of larger sizes. To this end, and differently from conventional techniques, according to the present invention the higher order statistical processing of the invention identifies coherence relashinoships between the results (coherent relationships/subnetworks) of lower order statistical processing conducted according to the invention, and either another lower order statistical processing conducted according to the invention (e.g. 1^st order pairwise coherence, or higher obtained in previous computational iteration) or one of the m signals which are to be modeled. Differently from 1^st order (pairwise) relations between pairs of signal, the higher order statistics can identify the relations/coherence between signals which are not coherent to 1^st order. For example, processing to order 0=2 may characterize connections (coherency) between two 1^st order networks by determining the coherent relationships between two pairwise coherent relationships obtained in the 1^st order and thereby yield a second order connectivity network of L=4 regions; or it may characterize the connections (coherency) between a 1^st order network and one of the signals itself (which may be considered here as a 0^th order network), and thereby yield a second order connectivity network of L=3 regions. Likewise, processing to order 0=3 may characterize coherent relationship (coherency) between two 2^st order networks (to yield 3^rd order network of size L=8), or between a 2^st order network and a 1^st order network (to yield 3^rd order network of size L=6), or between a 2^st order network and a 0^th order network (a signal), (to yield 3^rd order network of size L=5).

To this end, the order O of the network actually indicates the number of computational iterations required to devise the coherent relationships in the network. In this regards, it should be understood, that the computational iterations are not necessarily performed with actual repetition/iterations of a certain computation, but may be encapsulated in a single signal/data processing operation/computation, e.g. implementing a single "formula" presenting the results of such "iterations" to a desired order. Thus processing to order 0=2 may s characterize up to quadrupole networks (networks of m = 4 regions); processing to order 0=3 characterizes connections between a quadrupole networks and a network of same or lower order to and is therefore capable of characterizing up to octupole networks (networks of m = 8 regions). Accordingly, in general signal processing to order O carried out according to the technique of the present invention is capable of characterizing networks of up to size L = 2°.

Operation 110 of method 100 of the present invention includes the provision of n signal components S"(T) (namely n signals/time-series data) acquired/obtained from a plurality of n distinct monitored regions of the monitored entity (e.g. Regions Rl to R4 in Figs. ID and IE and/or regions Rl to R8 in Fig. IF). Here S"(T) denotes the signal/time-series S obtained from the n^th region R_n as a function of time T.

To this end system 200 includes a Data/Signal Receiver Module 210 which is configured and operable carrying out operation 110 of method 100 by connecting to a measurement module 282 and/or to a Measurements-Data/Signal Storage module 292 who stores such measurement signal. Thus, Data/Signal Receiver Module 210 provides the system 200 with data indicative of the n signal s/time-series(s) S"(T) acquired from a plurality of n distinct monitored regions {R} of the monitored entity E. The measurement module 282 may be for example an MRI, EEG, electrophysiology, and/or MEG system, and/or seismograph system and/or stock market data retriever module, and/or other measurement system, all depending on the discipline to which the technique of the present invention is applied.

In various implementations the measurement module 282 and/or Measurements- Data/Signal Storage module 292, being the data source of the signals/data S"(T) are configured as integral part of the system 200. In other implementations the Data/Signal Receiver Module 210 is configured for connecting/communication with such data source(s) (e.g. being local or remote) that do not form part of the system.

In optional operation 120 of method 100. at least a subset of m of the acquired n signal components S"(T) are processed to extract relevant signal portions thereof. The size of m (being the number of signals) may generally depends on the order O whereby the maximal size of the functional network/subnet that can be characterized by the processing is m~2°

Generally the technique of the invention is directed to determine connectivity networks/subnets of high order 0>2 between 4 or more regions of the monitored entity. Accordingly, the number n of the signals/regions, as well of the size of the subset m of signals/regions that are processed to determine the connectivity network between them is generally at least 4. This enables to determine high order connectivity networks which may serve as accurate markers/indicators for diagnosing phenomena/disorders in the monitored entity/system. In this regards it should be understood that the disorder/phenomena indication provided by high order connectivity network/model is substantially more reliable and accurate than low order models (e.g. of 0=1) which indicate only pairwise connections/couplings between regions/signals. To this end one or more subsets of m=2° (where 0>2 is the order of the connectivity network processing) are elected from the provided n signals (i.e. n>4 ; m>4; and n≥m) and respectively processed to obtain one or more respective subnets presenting a connectivity network model of the regions of monitored entity/system. In the following a detained description of the technique for determining a connectivity network model of each subnet of a subset of m>4 signals is described. This may be performed one or more times to obtain one or more subnet models for one or more subsets of m regions in the monitored entity.

In some embodiments of the invention, the elected subset(s) of m signals are further processed to extract relevant time and/or frequency portions thereof, which are relevant for the characterization of the phenomena/disorder of interest in the monitored entity E.

For instance, it may be a-priory known that for characterizing phenomena events, such as earthquakes, epileptic seizure, heart attack, certain predetermine time windows of the m signals, (e.g. preceding/proceeding or during the event) should be processed. Accordingly, in optional operation 122 time slot data indicative of a time slot of interest of the m signals is provided/obtained. The time slot data may be indicative of the timing of the time window/slot. In this regards the time slot data may be obtained by any of the following:

(iv) providing event timing data indicative of a timing of occurrence of event associated with the disorder/phenomena to be characterized and utilizing the event timing data to set the time window timing;

(v) setting the time window duration based on a type/time scale of said disorder/phenomena of interest- e.g. earthquakes/seismological phenomena should be analyzed in an order for tens of minutes/hours or even days; epileptic seizures in an order of tens of milliseconds; information transfer in the brain by MEG data in an order of tens of milliseconds, etc.;

The above timing and/or duration data may be a priory stored as a configuration parameters of the system/method for processing specific phenomena, or it may be provided as input data (e.g. from and operator of the system 200A) indicating of at least one of the: (a) time window timing, (b) time window duration, and (c) a time resolution of the signal to be processed.

Additionally or alternatively, it may be for instance a-priory known that in characterization of certain phenomena, certain predetermine frequency bands of the m signals are of more interest than others. For instance EEG data is expected at specific frequency bands generally known as delta, alpha, beta and gamma frequency bands. Accordingly, in optional operation 124 provided is frequency band data indicative of one or more frequency bands of interest in the m signals. The frequency band(s) may be a priory determined as configuration parameters of the system 200A, or they may be determined in accordance with the type of the phenomena/disorder that is to be characterized; and the relevant frequency band(s) may be provided as input data from the operator of the system.

In the following description generally the processing of the entire signals (e.g. including all its frequency and time portions is described. However it should be understood that generally the similar processing may be conducted several times for one or a plurality of frequency band and time window corresponding portions of the subset of m signals. Per each time window and frequency band, the connectivity network may be characterized, and the plurality for connectivity networks for the different frequency bands and time windows may be further processed by the system, or considered by an operator of the system, in order to assess the phenomena of interest.

Thus, in optional operation 126 the subset of m signals are processed to extract respective frequency-time domain portions of the signal corresponding to the provided time slot(s) and frequency band(s) of interest. ;(optionally utilizing time-frequency wavelet transform). In the following described in details is a specific implementation of frequency-time domain processing carried out utilizing the so called wavelet transform. However, it should be understood that the invention is not necessarily limited to use of wavelet transform, and that possibly other frequency-time domain processing such as Fourier transform, possibly carried out on time windows of the m signals, may be used for extracting relevant time window and frequency band portions of the signals. Thus, in optional operation 126, the m signals may be transformed via a proper time-frequency transformation, e.g. wavelet transform, so that m respective transformed representations (e.g. complex/wavelet representations) of the m signals are obtained.

information of similar frequency.

in Eq. (1) below. Also the selected wavelets and the form of the wavelet representation are exemplified in more details with reference to in Eqs. (2) and (3) below.

It should be noted that the term wavelet representation used in the following should be understood as general designating any type of frequency-time domain portions of the m signals (e.g. obtained in operation 120 above), which may be obtained by utilizing any suitable technique (not necessarily wavelet transform), and which may include specific time-frequency portions of the m signals or alternatively include the of the entire m signals of the complete information thereof.

To this end, system 200A includes an optional Frequency-Time domain signal processor Module 220 adapted for extracting the relevant time window and frequency band portions from the m signals. The Frequency-Time domain signal processor Module 220 may for example include a TimeSlot-FrequencyBand provider module 223 configured and operable for carrying out operations 122 and/or 123 of method 100 described above, for providing data indicative of the relevant time slots and/or frequency bands to be extracted from the m signals. Also optionally the Frequency- Time domain signal processor Module 220 may include a Wavelet Transform Module 226, which is configured and operable carrying out operation 126 of method 100 as described above. The Wavelet transform Module 226 may be a computerized module/system, implemented by software and/or hardware, and may be connectable to, or may include processing unit (e.g. a central processor (CPU) and/or a digital signal processor (DSP)) and to a memory) operable for performing the operation 120.

Operation 130 includes sub-operations 132 and 134 which allow to second order O = 2 or higher statistical association between at least one subset of m regions/signals (e.g. between m portions or wavelet representations of the signals), and optional operation 136 enabling to obtain statistical associations of even higher orders.

In operation 132 the m signals (or portions thereof) are processed to determine low order (0=1) functional connectivities between them. This yields pairwise connectivities between pairs of the m regions/signals. The low order processing is performed by determining respective pairwise coherences between the signals of respective pairs of the m signals. This yields data indicative of the low order functional connectivities between the respective pairs of the m regions of the entity/subject being monitored.

In turn, up to a 1 ^st statistical order, regions for which exist a pairwise coherence between the functional representations of their respective signals, are considered/classified as functionally connected regions. Otherwise, to the 1 ^st statistical order computed in operation 132, such regions may be classified as disconnected.

Operation 134 of method 100 is carried out in order to identify/reveal functional connectivity of higher statistical order(s) (e.g. 2^nd and/or higher) between the m regions of the subset. This is achieved in operation 134, by processing the functional connectivities (e.g. FC[l-2] and FC[3-4]) of the lower order, which are obtained in preceding operations of method 100 (e.g. in operation 132, or in preceding iterations of operation 134 - in case more than one iteration is conducted to obtain connectivity networks of order higher than 2), to determine coherent relationships between the functional connectivities of the lower order (e.g. identifying coherent relationships between FC[l-2] and FC[3-4] or identifying coherent relationships between FC[l-2] and R4 as shown in Fig. IE).

To this end, according to the present invention, existence of functional connectivity (e.g. existence of a coherent-relationship) between a pair of the lower order functional connectivities (e.g. of order O) that are obtained in the preceding operation (e.g. 132) is indicative of functional connectivity of a higher order (e.g. of order 0+1) between a subnets of a plurality of L = 2⁽⁰⁺¹⁾ regions of the m subset regions (whereby the identity of the L regions is that of the regions participating in the pair of low order functional connectivities from which the higher order connectivity is obtained). Accordingly, operation 134 provides for obtaining a high order connectivity network model of the m regions which may serve as a marker/indicator of a disorder/phenomena expressed in the monitored entity/system. Advantageously, operation 134 allows to efficiently model the connectivity network between large number of regions, by using statistical processing of 2^nd or higher statistical orders. Operation 134 is described and exemplified in more details below with reference to Eqs. (7) and (8).

Optionally, as illustrated in operation 136 in Fig. 1A, in some embodiments of the present invention operation 134 is repeated/iterated for one or more additional times to obtain functional connectivity of even higher orders (e.g. 3^rd order - 8 region subnetworks). In other words, each iteration multiplies (doubles) the statistical order, for example the second iteration order lead to determining connectivity subnet of 4 signals/time series, (2^nd order statistics). Indeed, by carrying out operation 134 to determine coherent relationships between the functional connectivities of order O (e.g. 0=2) obtained in a preceding iterations of operation 134 a higher order coherent- relationship between the monitored regions (i.e. of order 0+1) can be reveal ed/i dentifi ed .

Thus, according to some implementations, when the subset size m is 4, the technique for determining complete functional connectivity network for the subset may terminate after operation 134 (this will yield single subnet of size L = m = 4). However, in case a determining complete functional connectivity network for subsets of m>4, the method 100A may further include 136 for repeating operation 134 for one or more additional 0-2 times in order to obtains a functional connectivity network of further higher order of O between the subset m regions, where O is an integer number greater or equal to log₂(m). In case a complete functional connectivity network model for all the m regions is needed, then operation 134 may be repeated for 0-1, to yield a single subnet of size L = m of up to 2° regions.

For example, Fig. IE shows a model of two functional connectivity networks of order 0=2 between the regions Rl to R4, similar to that of Fig. ID only that here operation the operation 134 was further carried out once. In the left example shown in the figure operation 134 had revealed, a functional connectivity of a second order FC2 [[l-2]-[3-4]] linking the first order functional connectivities FC[l-2] and FC[3-4] between the respective regions R1-R2 and R3-R4. In the right network example shown in the figure operation 134 had revealed, a functional connectivity of a second order FC2 [[1-2J-4] linking the first order functional connectivities FC[l-2] and the region 4. To this end, operation 134 reveals that their exists an actual functional connectivity between the regions Rl and R4, which are considered, disconnected to the first order (see Fig. ID)

Fig. IF shows to examples of models of the functional connectivity network between eight regions Rl to R8, which is revealed by method 100 by carrying out operation 134 twice (i.e. iterating/repeating twice the operation 134 as indicated in the optional operation 136 of the method). In the left network exemplified in the figure, not only that second order functional connectivity, FC2 [[l-2]-[3-4]] and FC2 [[5-6]-[7-8]] are revealed, but also a 3^rd order functional connectivity FC3 is revealed/identified which links the 2^nd order functional connectivities. In the right network exemplified in the figure, second order functional connectivity, FC2 [[l-2]-[3-4]] is revealed and also 3^rd order functional connectivity FC3 is revealed/identified which links the second order functional connectivity FC2 [[l-2]-[3-4]] with a first order functional connectivity FC[6-8]. In both examples a model of the connectivity network between the regions Rl to R8 is obtained .

In this regards it should be noted that the technique of the present invention is advantageous over conventional techniques for modeling connectivity networks, such as nervous system's connectivity, since the conventional techniques are based on the statistical identification of pairwise interactions between respective pairs of regions, and in this sense can only identify connectivity between regions which present statistical association/dependence of a first order between their signals. For example utilizing the conventional techniques regions Rl and R4 in Fig. 1C would not be considered functionally or directively connected, unless exists an additional 1^st order functional connectivity, such as FC[2-3], which provides for linking the regions Rl and R4 with first order functional connectivities connecting the immediate regions between them. In other words, according to the conventional techniques, in order to identify the functional, or directed connectivity between regions (e.g. Rl and R4) of a nervous system or other entity, there should exist a complete path of pairwise (1^st order) statistical associations/correlations linking these regions Rl and R4 (e.g. possibly via one or more intermediate regions such as R2 and R3 between them).

Furthermore, it should be noted that technique of the present invention is advantageous over conventional techniques since is allows to identify high order statistical couplings between pluralities (>2) of simultaneous signals portions (e.g. associate with the same time window and/or with the same frequencies) this reveals simultaneous interactions between the plurality of m>2 signals/time-series).. Another advantage of the technique of the present invention, is that by carrying the analysis on different time windows and/or frequency bands while controlling/adjusting the parameters of the time windows or frequency bands, various additional characteristics of the phenomena/disorder can be determined (such as the duration of the phenomena, and/or its prominent frequencies).

This presents a deficiency of the conventional technique since the regions, e.g., Rl and R4, may be still in practice functionally or directively connected, even in cases where a complete path of pairwise (1^st order) statistical associations/correlations does not exist or is not identified between these regions (e.g. between Rl and R4) the regions may still be functionally connected (e.g. possibly via other regions, such as regions R9 and RIO in the figure, which may not be monitored and/or from which signals are not acquired). To this end, the technique of the present invention allows to identify such directed and/or functional connection between regions of the nervous system which would be considered disconnected according to the conventional techniques.

Overcoming this deficiency of the conventional techniques is achieved according to the present invention by using the higher order statistical analysis of the signals acquired from the monitored regions. This allows identification of statistical dependence/association between the activities of sub-networks in the connectivity model which would be considered otherwise disconnected. This has tremendous significance in the ability to recognize and diagnose various phenomena/disorders which might not be otherwise identified/revealed in the subject/entity that is being monitored.

Thus operation 130 yields a model of a functional network indicative of the functional connectivities between a plurality of more than two monitored regions (n regions) in the subject/entity being examined. For instance, carrying out operation 130 of method 100 may reveal models of the functional connectivity network, such as those illustrated in Figs. IE and IF, between plurality of regions (4, 8 or even more regions) of the examined subject/entity.

Typically, relatively efficient processing of the signals from the n regions to determine their functional network is achieved in operation 130 when the number m of the subset of regions to be included in the functional network is an exponent of 2 by the integer number O being the order of the processing (O also represents the number of repetitions of sub operation 134 plus one). In other words, relatively efficient processing may be achieved when the number of nodes/regions m = 2°. This is because each of the O-l repetitions of the sub-operation 134 doubles the number of regions that can participate in the resulting model of the connectivity network. To this end, when sub operation 134 is repeated (actually and/or merely mathematically) once (2^nd order), and the number n of regions is 4 regions. Namely, in this case operation 130 results with a 4-region functional connectivity network. In the same way, when sub operation 134 is repeated twice, the statistical order O is 3), m = 8, so in this case operation 130 results with an 8-region functional connectivity network. In the same manner larger connectivity networks may be efficiently modeled by utilizing processing of higher statistical orders.

With regards to operation 130, it should be understood division of operation 130 in the description above into the sub-operations 132, 134 and 136 is made only in order to clearly describe the method of the present invention. Indeed, such division makes understandable the rational underlying the functional modeling technique implemented in operation 130, however it should be understood that practically the entire operation 130, and optionally both operations 120 and 130, may be implemented simultaneously within a single signal processing operation. As will be readily appreciated by those versed in the art, and as is also clear from the example of Eq. 7 below (which encapsulates the entire operation 130 for a particular example of a 4-region functional connectivity network), the entire operation 130 and optionally also together with preceding operation 120, may be implemented together, by a single operator/module implementing in a single signal processing operation. In the same way, modeling of higher order networks (e.g. of 8 and/or 16 nodes) can also be implemented.

System 200A includes Functional Connectivity Mapping Module 230, which is configured and operable carrying out operation 130 of method 100A as described above. Functional Connectivity Mapping Module 230 may be a computerized module/system, implemented by software and/or hardware, and may be connectable to, or may include processing unit (e.g. a central processor (CPU) and/or a digital signal processor (DSP)) and to a memory) operable for performing the operation 130.

Turing now to Fig. IB, a flow chart of a method 100B for determining an hierarchical and/or directed connectivity network between the subset of m regions is illustrated. According to some embodiments the method 100B, utilizes the output functional connectivities obtained by the technique of the method of 100A (or possibly in case equivalent technique for determining high order connectivity networks will be found, utilizing the results of such equivalent technique), and further processes the functional connectivities (e.g. processing the their characterizing coherence relationships) to determine the directive connectivity network between the m regions (e.g. the directed connectivity and/or the hierarchy between the connected regions).

Indeed, for some applications of the technique of the present invention, obtaining a model/map of the functional connectivity network between regions, may not be sufficient to identify/diagnose the existence of a phenomena/disorder in the examined subject. More specifically, for some applications accurate diagnostics/identification of a phenomena requires a map/model of hierarchical connectivity network between the examined regions (e.g. such as that illustrated/exemplified in Fig. 1G) and/or even a model of the directed connectivity network (e.g. such as that illustrated/exemplified in Fig. 1H). In this regards, it is noted that a model of the functional connectivity network merely indicates which regions are functionally connected/coupled (e.g. being activated in correlation with one another). A model hierarchical connectivity network, such as that of Fig. 1G, indicates, in addition to the functional connectivity, also the order/hierarchy of the activation of the regions, while not necessarily indicates the direction/causality of the activation. A model of the directed connectivity network, such as that of Fig. 1H, indicates, in addition to the hierarchical connectivity, also the direction/causality of the activation of the regions.

Therefore, optionally, according to some embodiments of the present invention, the method 100B for modeling the connectivity network further includes operation 140. Operation 140 provides for determining the temporal hierarchy (order) between the regions/nodes of an functional connectivity network model (e.g. a high order model such as that obtained by method 100A), and optionally also provides for determining the temporal directionality (the causality) of the modeled network. The latter thereby provides for determining/modeling the directed network between the regions.

In this regards, it is noted that a feature of the technique of the present invention is that it allows to determine the order/hierarchy (not necessarily the directionality) between the regions included in the modeled connectivity network. The order/hierarchy may be determined based on the phase relations between the subset of m signals S"(T) for which the functional connectivity network is determined. This is exemplified in more details below with reference to Eq. (9). More specifically, as shown for example in Eq. (9) below, according to some embodiments of the present invention, the phase relations between the m signals S"(T) (or portions thereof) are determined accordance with the relative phases of the functional connectivity's (coherences of the signals) determined in the functional connectivity network.

To this end, according to some embodiments of the present invention operation

140 includes processing/computing/determining phase relations between n-1 signals acquired from n-1 of the regions relative to the phase of the signal of another region (e.g. an η^Λ region) of the network, and determining the temporal hierarchy/order in accordance with said relative phases. Accordingly an hierarchical connectivity network indicative of the ordered relation between the regions of the network is determined. For instance, Fig. 1G illustrates the hierarchies hi, h2, h3 and h4 assigned to the respective regions R2, Rl, R4 and R3 in this way. This means that the activation order of the regions may be either from R2 to R3 as follows: R2→ Rl→ R4→ R3 , or from R3 to R2 as follows R3→ R4→ Rl→ R2.

Optionally, in some implementations, operation 140 further continues to determine the directed connectivity network between the regions. In such implementations operation 140 further includes providing an assumption/reference data which is indicative of the temporal directionality/causality of at least two of the regions in the network. Such assumption/reference data may be obtained for example from a from a suitable data storage storing a reference model/data about the type of the subject/entity which is being examined. Such assumption/reference data is used together with the hierarchical connectivity network (which may be determined based on the phase relations as indicated above), in order to determine the complete directional/casual hierarchy/order between the regions of the network and thereby model the directed connectivity network which is formed by the monitored regions. Fig. 1H illustrates an directed connectivity network modeled in this way whereby the arrows indicate the causality relation between the regions.

Accordingly, optionally, system 200A may further include an Directed

Connectivity Mapping Module 240 configured and operable for carrying out operation 140 of method 100 as described above, to determine the hierarchal and/or directed connectivity network model of the monitored entity/system E. The Directed Connectivity Mapping Module 240 may be a computerized module/system, implemented by software and/or hardware, and may be connectable to, or may include, a processing unit (e.g. a central processor (CPU) and/or a digital signal processor (DSP)) and to a memory) operable for performing the operation 140.

Optionally, according to some embodiments of the present invention, the method 100 further includes operation 150, which is carried out to determine the type of the connectivity network between the monitored regions, and/or more specifically to classify the network a linear network, non-linear network and/or a disconnected network.

In some implementations the linearity of the network is determined based on a linear relationship (see for example Eq. (12) below) between the relative phases of pairwise coherence (Eq. (4) below) of respective pairs of the n wavelet representations and the phase relations (Eq. (9)) between m-1 of the regions relative to the m"¹ region.

For example, a connectivity network composed of m = 4 regions can be constructed in accordance with method 100A above, to obtain after operation 130, a 4- region- functional connectivity network (or in more general m-region- functional connectivity network if higher order statistics is employed in operation 130). Alternatively or additionally, another technique to compose a connectivity network of m = 4 regions, is by using the six possible pairwise combinations of the 4 regions to compute the pairwise coherences (e.g. as shown in Eq. (4) below) of the signals of these six combination pairs. The latter technique yields a so called "Joint-pairwise- coherences" of the four regions. Thus according to some embodiments both the technique for constructing the connectivity network are carried out to determine both types of connectivity network models, the m-region-functional connectivity network and the Joint-pairwise-coherences composed of m (e.g. m=4) regions. Then both types of connectivity networks models are compared in order to determine the type of the connectivity network of the entity/system E. In case that both analyses: the m-region- FC model of the present invention (operation 130 above), and the joint-pairwise- coherences model result with a significant functional connectivity (FC), it implies that the connectivity network is a dynamically stable network - which is referred hereafter as a linear network. In only the m-region-FC model yields significant functional connectivity, it implies that the network is a combined network. In case only the joint- pairwise-coherences model results with a significant functional connectivity, it implies that significant coherences between different pairs of regions occur at different times, suggesting that the entire network does not work together or that the regions are functionally connected through higher coherences with higher (>4) number of regions. Such network is referred hereafter as a disconnected network (note that the disconnected definition here differs from the definition used in graph theory).

To this end, optionally, system 200A further includes a Network Type Classifier Module 250, which is may be configured and operable carrying out operation 150 of method 100 as described above, in order to determine the type of the modeled network. Network Type Classifier Module 250 may be a computerized module/system, implemented by software and/or hardware, and may be connectable to, or may include, a processing unit (e.g. a central processor (CPU) and/or a digital signal processor (DSP)) and to a memory) operable for performing the operation 150.

Typically the functional and/or hierarchical and/or directed connectivity network models and/or the type of the network, which are obtained in operation 130 and possibly 140 and/or 150 above are stored in a data storage for further use, and possibly stored in association with the subject for which these models were constructed. To this end, system 200A may further include, or be associated with, an Entity's Connectivity Model Data Storage 272 configured and operable for storing the such connectivity network model data for the subject. According to some embodiments, the functional and/or hierarchical and/or directed connectivity networks and/or the type of the network, which are obtained in operation 130 and possibly 140 and/or 150 above, are further used to identify/diagnose a disorder/phenomena/condition expressed in the examined subject/entity. To achieve that, optionally method 100B includes operation 160 which includes: (i) providing reference data which may for example include one or more reference models of the connectivity network between corresponding regions Rl to Rn in normal (e.g. healthy) and/or in abnormal (e.g. disordered) entities (subjects/systems); and (ii) comparing the connectivity network model (e.g. the functional and/or hierarchical and/or directed connectivity network) obtained for the examined subject in operation 130 and/or 140 with the reference model of a normal and/or abnormal entities (subjects/systems), and (iii) determining a similarity between the connectivity network model of the examined subject and the reference model to assess whether the phenomena/disorder which is looked for is expressed in the entity being examined or not (namely classifying the entity/system to be normal or abnormal.

For instance Fig. II schematically illustrates in self-explanatory manner an example of a reference data indicative of a connectivity network reference model of normal entities/subjects expressing no particular disorder/condition/phenomena. Figs. II and Figs. 1 J are self-explanatory illustrations of two examples of connectivity network reference models of entities/ subjects expressing/having different types of di sorders/ conditions/phenomena.

The reference data (e.g. the reference connectivity network models which may include for instance: (i) a plurality of connectivity network models of a plurality of subjects expressing no particular disorder/conditi on/phenomena, and/or (ii) an average/nominal/representative model of such "normal" subjects, and/or (iii) a plurality of connectivity network models of a plurality of subjects expressing certain particular disorder/condition/phenomena; and/or an average/nominal/representative model of such "abnormal" subjects, may be stored in a reference data repository, for example the Reference Data Storage 274 which may be included and/or associated/connectable with the system 200A.

Optionally, system 200A includes Phenomena/Disorder/Condition Identifier module 260 (e.g. including a Connectivity Model Comparator). The Identifier module 260 may be adapted to identify/diagnose a phenomena/condition of the subject by comparing the connectivity network model of the subject (as obtained in any one or more of operations 130 to 150 and possibly stored in the data storage 272) with the reference connectivity models (e.g. stored in the Reference Data Storage 274 which stores Phenomena associated Reference Connectivity Models) to determine/diagnose whether the entity/subject E expresses a certain phenomenon/condition/disorder of interest or not. This may be achieved carrying out operation 160 of method 100B for determining the similarity(ies)/difference(s) between the subjects' own connectivity network model and connectivity network models of abnormal subject and/or of normal subjects expressing/not-expressing the disorder.

The Phenomena/Disorder Identifier module 260, which is configured and operable carrying out operation 160 of method 100B, may be a computerized module/system, implemented by software and/or hardware, and may be connectable to, or may include processing unit (e.g. a central processor (CPU) and/or a digital signal processor (DSP)) and to a memory) operable for performing the operation 160.

Optionally system 200 further includes a Reference Data Generator 275 that is adapted to receive external/independent data about the condition of the subject E and utilize this external data and the connectivity network model determined for the subject E in any one or more of operations 130 to 150 to update the reference data stored in the Reference Data Storage 274.

In the following provided are several specific exampled for the application of the technique (method and system) of the present invention in the field of neuroscience. In this field the technique of the present invention is particularly advantageous since the methods 100 and 100B allow diagnosing nervous system disorders/phenomena/conditions in a noninvasive manner. The technique of the present invention provides a novel technique by processing and/or transforming signals measured from the individual's nervous system, to obtain a connectivity network model of the nervous system that allows identification of clinical markers indicative of one or more disorders or types thereof.

The obtained model may be characterizing of the connectivity network of the individual's nervous system and more specifically characterizing the connectivity between plurality of regions of the nervous system, in between the functional and according to some embodiment the directed connectivity (network) of the individual's nervous system, which for that matter may refer to the central nervous system of the individual (e.g. its brain) and/or his peripheral nervous system.

According to the invention the modeling of the nervous system is achieved by processing signal components, such as MRI (e.g. fMRI), EEG, MEG, and/or electrophysiology, acquired from a plurality, n, of more than two distinct regions of the individual's nervous system, and determining "simultaneous" network connectivity model between the signal components from these plurality of m > 2 regions.

Below provided is a detailed mathematical description of the method 100A and 100B of the present invention implemented particularly for characterization of a 4-region connectivity network. As will be readily appreciated by those versed in the art, generalizing the mathematical description to higher orders (larger networks is readily possible by implementing the operation 136 described above, on the description of the 4 regions network connectivity models, which is provided in the following in mathematical terms.

For clarity, the non-limiting examples provided below refers particularly to the processing of nervous system signals which are obtained via resting-state functional MRI. However it is understood that the technique of the invention as described above and below is not limited to resting-state functional MRI and can be implemented on other coupled time-series data in varies disciplines such as seismology, or economy. In neuroscience, it can be used with stimulus-driven fMRI, electrophysiology, Magnetoencephalography (MEG) or Electroencephalogram (EEG).

BOLD fMRI time-series are generally ergodic (in a weak sense). Therefore, time frequency analyses and spectral coherence are appropriate ways to summarize functional connectivity in a compact and efficient fashion. A pair of time-series is considered coherent if they have a constant relative phase for a given time and frequency windows. In the following examples it is shown that by using wavelet analysis and calculating connectivity according to spectral coherences, connectivity network models of multiple time-series (each corresponds to a different region) can constructed/generated and provide their directed functional information.

For a given frequency window, the network's coherence is defined by phase- relations of all time-series in the network, namely, a combined coherence. This is in contrast to networks defined by all possible pairwise coherences among network time- series, namely, a linear coherence. As indicated above, the technique of the present invention analysis described below with reference to BOLD signals, can be applied to electrophysiological and other imaging modalities that are not confounded by variable hemodynamic delays. Consequently, assuming that time-lags of BOLD signals correspond to latencies of information flow, BOLD signal phases are used to infer the temporal hierarchy with no need for other assumptions. Directed hierarchy (corresponding to directed functional connectivity) is obtained here by assuming a simple relation (e.g., region A, out of the four, precedes region B), with no need for further statistical models.

highest frequency defined by the band-pass filter) and t the time (shift).

The complex wavelet coefficients can be written in terms of their amplitude and phase:

As indicated above Eqs. (1) to (3) present an example of the transformation/processing operations which are applied to the time series signals data in operation 120 of method 100, and/or by module 220 of the system 200 of the present invention. Mathematically, a Pairwise fiinctional connectivity (pairwise FC) may be represented as follows:

and phases, respectively, of the BOLD signal wavelet coefficients in Eq. (3).

As indicated above Eq. (4) present an example of the transformation/processing operations which are applied to the time series signals data in operation 130 (specifically 132) of method 100, and/or by module 230 of the system 200 of the present invention.

Since resting-state data is not time-locked to any particular event, averaging over time may optionally be performed at the last stage to yield a FC that is dependent only on frequency and not time:

To obtain a single (complex) FC value, averaging over frequencies can also be done, yielding a FC value which is similar to our previous definition using Fourier space analysis (see Ref [10]), and whose real part corresponds to the Pearson' s correlation coefficient:

values and Pearson's correlation coefficients across all pairs was 0.9577 (p value practically zero). On average, the values obtained by Equation 6 were higher than Pearson's correlation coefficients by 27%.

Since Equation 3 for the BOLD signal and Equation 4 for pairwise FC have the same form, in Equation 4 can be replaced by the FC of two regions

independent ways to do these insertions. The 4-region-FC expression is defined as the multiplication of these three possible combinations:

Fig. 2 shows an example of 4-region-FC for four subcortical regions that exhibit FC at an intermediate frequency range. More specifically, Fig. 2 shows an example of the technique of the present invention implemented fort to obtain the nonlinear interaction between 4 time-series (named 4-region-FC). Sections (A) to (D) Wavelet amplitudes shown in time-frequency space, with time in the X-axes and frequency in the Y -axes, for 4 time-series signals. The signals were obtained by functional MRI acquired during resting state from 4 different regions in the human brain of a single volunteer. Sections (E) and (F) show respectively the amplitude and phase of the time- frequency representation of the combined network composed of the 4 signals. It points to time-frequency windows in which the combined amplitude is high and the combined phase is constant. The combined amplitude and phase are given in equation 7. In (E) on the right side, the average along time pf the real (solid white) and imaginary (solid red) part of the wavelet representation of the whole network is shown as a function of frequency. As indicated above Eq. (7) present examples of the transformation/processing operations which are applied to the time series signals data in operation 130 (specifically 134) of method 100, and/or by module 230 of the system 200 of the present invention.

with N = k + I. For example, the left network shown in Fig. IE is formed with k and / being 2 and 2 respectively, while on the right network in this figure, k and 1 are 2 and 1 respectively. Both networks are of order 2. Similarly on the left network shown in Fig. IF is formed with k and / being 4 and 4 respectively, while on the right network in this figure, k and 1 are 4 and 2 respectively. Both networks are of order 3. Note that k does not have to equal to / and either can be odd or even numbers thus enable N to be any number. As indicated above, according to some implementations of the present invention the phases of differences between coupled time-series functions are used to obtain directionality of the modeled connectivity network.

For 4-region-FC of Equation 7, the three independent phase expressions of Equation 8 enable to obtain the phases of three BOLD signals (any three) as a function of the forth:

As indicated above, in some implementations of the present invention operation 150 of method 100 is implemented to process the signals are further processed/transformed in order to determine the type of the connectivity network between the regions. Namely to classify the network for example to one of the following types: Linear, combined and disconnected network_.

10. The phases for a network calculated from six pairwise coherences are:

stable network. Such a network is referred hereafter as a linear network. If only the FC obtained by the 4-region-FC calculation is significant, the phase-relations given in Equation 7 must hold but Equation 11 is not satisfied. This type of a network is referred hereafter as a combined network. If only the joint-pairwise-coherences result with a significant FC (Equation 11), it implies that significant coherences between different pairs of regions occur at different times, suggesting that the entire network does not work together or that the regions are functionally connected through higher coherences with higher (>4) number of regions. Such network is referred hereafter as a disconnected network. Note that this definition differs from the definition used in graph theory for fragmentation.

To this end, Eqs. (10,11) presents an example of a technique for classifying networks, which may be implemented in operation 150 of method 100, and/or by module 250 of the system 200 of the present invention.

It should be noted that the Example for 4-region coherence connectivity network described above, including its dependency on coupling strength and noise level, was verified by the inventors by utilizing computer simulations (see Figs. 7 and 8 described below). Also, the validity of this technique was further demonstrated based on resting- state functional MRI data of 40 young healthy subjects to characterize the ventral visual system, motor system (see for example Figs. 9 to 12 described below).

Reference is now made together to Figs. 3A and 3B which show flowchart 100C and block diagram 200B schematically illustrating a method and a system according to an embodiment to the present invention for identifying /determining a region of source of a disorder/phenomena expressed by a certain entity (system/subject). It should be noted that common reference numerals presented in these figures and in Figs. 1A to 1C, which are described in detail above, pertain to like/similar elements/modules/operations, as described above with reference to Figs. 1A to 1C. Accordingly, for clarity and conciseness, description of these elements/modules/operations need not be repeated in details in the following.

As described above with reference to Figs. 1A to 1C, the present invention provides a novel technique for processing signals (sequential series) obtained from a plurality of regions of a monitored entity/system processing these signals to modeling these respective regions as a network having functional connectivity, and possibly also hierarchical and/or directed connectivity, between the regions. This is achieved by utilizing high order (0>1) statistical tools to infer none-linear coherence between the subset(s) of plurality of m regions from which m respective signals are obtained. To this end, the functional connectivity between regions is determined by determining the high order coherent relationships between the signals acquired from these regions respectively, while hierarchy between the regions, is determined utilizing the phase relations between the coherent relationships determined for the functional connectivity. The connectivity model/network obtained in this way, may be used to identify/diagnose various phenomena/disorders in various types of monitored entities.

However, in some cases, for example when attempting to diagnose and/or treat disorders/phenomena, that are associated with seizure like behavior, such epilepsy/, earthquakes, or other, there arises a need to identify not only the connectivity network behavior of the regions of the monitored entity, but also determine/identify the regions that are most likely the source of the disorder (the source of the seizure like behavior).

To this end, method 100C and system 200B, provides a novel technique according to an embodiment of the present invention for identifying/determining region(s) of the monitored entity, which are most likely a source of the disorder or seizure associated with it. In some cases (for example when considering epileptic disorders), identification of such source regions, facilitates the efficient and accurate treatment of the disorder.

The proposed technique exploits the ability of the method 100A and 140 described above for inferring temporal hierarchy in networks of subsets of m coupled regions/signals in order to determine which of the n regions from which the n signals are obtained is likely the source of the disorder.

The method 100C includes operations 170, in which out of the n signals (regions) obtained from the monitored entity, K groups of m-1 signals, are selected (whereby per each group one signal of the m-1 signals is set as the reference signal). As will be clarified below, preferably according to some embodiments of the present invention the selected groups of m-1 signals and their associated selected reference signal are remained fixed in the further processing (at least per frequency band and time window for which the processing is conducted). Based on the selected K groups of m- 1, K*n subsets of m signals each are constructed, by adding each of the n signals to each of the K groups, such that per each of the n signals there are formed k subsets of m signals each, which are constructed with the different groups of m-1 signals.

Operations 180 and 190 includes carrying out operations 100A and 140 described above, per each one of the n*K subsets of m signals each. In this way hierarchical connectivity network model is obtained per each one of the K*n subsets. As indicated above, the hierarchical connectivity model is indicative of the relative phases between the m signals of each subset. These, as will be explained below play significant role in identifying the most likely regions of the n regions that are candidates to be the source of the disorder. This is because the phases of the signals sourcing the disorder, generally precede the phases of all other signals (the come temporally earlier). However, indeed, due to the cyclic properties of phase, it is difficulty, and practically impossible to identify with certainty, from a single network connectivity model, the phase of which of the n signals precedes all others, (this is also in part due to the large numbers of regions n~100s) which makes it difficult (computationally) to obtain a single complete model for such large network. To overcome this problem, the inventors of the present invention have devised a statistical method for inferring the candidate regions of the n regions, whose signals/phases precede most/all others of the n regions. This is based on the understanding that the phases of the regions, which are likely to be the sources of the disorder, precede the phases relative to other regions and therefore may accumulate (when averaged over a plurality of subsets/subnets) to a high value. On the contrary, the phases of signals from regions that are not the source of the disorder, may precede the phases of some reference (for some subnets/subsets), while proceed the phases of reference signals of other subsets, and thus on average be suppressed to nearly zero.

Accordingly, the operations 180 and 190 of method 100C includes the following:

- selecting one of the n signals to be a candidate signal associated with a source region of the disorder;

- utilizing the candidate signal with each one of the a plurality of k different groups of m-1 signals to construct a plurality of k subsets of m signals, each including the candidate signal;

- carrying out operations 100 and 140 on each of the plurality of k subsets which include the candidate signal to determine relative phases between the candidate signal and the reference signals of the respective subsets/groups (e.g. utilizing the phases obtained by Eq. 10 above); and

- summing the relative phases to obtain an average relative phase of the candidate signal. This may be repeated per each one of the n signals selected to be the candidate signal, and the candidate signal having the highest average relative phase may be selected as the most likely candidate for being the source region of the disorder. As indicated above, preferably, the operation is repeated for while forming the k subsets per each of the n candidates of the n signals, using the similar k different groups of m-1 signals each similar group is characterized by the identity of the its (m-1) signals of the subset and the identity of the reference signal therein. This yields suppression of average relative phases of candidate signals which are not associated with the source region of the disorder.

To this end, according to some implementations the method 100 may be more efficiently implemented as follows:

In operation 170 a group of m-1 signals is selected (e.g. randomly) from the total n signals whereby one of said m-1 signals of the group is selected as a reference signal. Then, operation 180 includes carrying out the method 100A and operation 140 discussed above for n times for n different subsets whereby each time a different one of the n signals is selected as a candidate signal and is joint to the group of m-1 signals to form together a subset of m signals on which the method 100A and operation 140 are performed. This yields n hierarchical connectivity networks of size m from which the relative phases of each of the n signals relative to the reference signal of the group of m- 1 signals can be determined (e.g. see Eq. 10 above).

Operation 190 includes repetition of the operations 170 and 180 above for k times (e.g. k may be for example of the order 1000) so that , k relative phases are determined per each signal of the n signals (being candidate signal) whereby the k relative phases are between the candidate signal and the k reference signals of the respective k groups of m-1 signals each.

Then the k relative phases of each candidate signal of the n signals are summed together to obtain for each candidate signal, an average relative phase.

As indicated above, it is expected that the average relative phase of one(s) of the n signals which precede most of the other n signals, will accumulate to substantial non zero value, while the average relative phases of others of the n signals will diminish/suppress.

Accordingly operation 195 identifying the most likely candidate signals to be the source of the disorder, those/that whose the average relative phase(s) obtained in operation 190 is/are the highest amongst the average relative phases computed for the n signals.

Turning now to system 200B, the operation 170 for selecting the K groups of m- 1 signals and constructing therefrom the n*K subsets of m signals each may be performed be implemented by the subset generator module 280, which may be a computerized module configured for carrying out this operation. The n*K subsets are then fed to sub-system 200A' which may be for example similar to system 200A illustrated in Fig. 1C, to determine hierarchical connectivity network model per each of the K*n subsets. The sub-system 200A' may be integral part of system 200B or it may be a separate system in communication with the system 200B. The hierarchical connectivity network models of the K*n subsets are indicative of the relative phases between each of the n signals and the references signals of each of the K groups. Accordingly, the Relative Phase Averager module may process the hierarchical connectivity network models of the K*n subsets to determine the average relative phase of each of the n signals (averaged/summed as described above relative to the K reference signals of the K groups). Finally, the candidate source region selector module 295 is configured for selecting and/or outputting the one or more of the n signals which have the highest average relative phase (e.g. exceeding a certain static or dynamic threshold) as the most likely candidates to be the source of the disorder.

Turning back to method 100C, it should be noted that according to certain embodiments the method is conducted for several J times for a plurality of frequency bands and/or time window portions of the n signals. In this regards, as indicated by operation 120 of method 100A described above, a connectivity network model may be determined per frequency band and/or per time window portions and/or per frequency- time portions of the signals. The inventors of the present invention have found that in certain seizure causing disorders, the seizure is expressed in several frequency bands. Therefore, it is preferable to conduct the method 100C described above, for several relevant frequency bands, and in this manner determine the most likely candidate region to consider as the source of the disorder/seizure, with improved accuracy and reliability. Accordingly, in some embodiments of the present invention, method 100C further includes operation 185 by which the above described operations 180 and 190 and optionally also 170, are repeated for several J times each with portions of the n signals which are associated with a different frequency band portion and/or a different time- window. Accordingly, in this case K*J relative phase are obtained per each of the n signals/regions (e.g. for the different K groups and the different J frequency bands). In this case, the averaging of the relative phases of each of the n signals in operation 195 is specifically conducted on the absolute value of the relative phases (this is because the sign of the relative phases is frequency dependent). Accordingly, as a result per each of the n signals there is obtained an average relative phase determined from over a plurality of frequency bands - thereby yielding more reliable indication as to the most likely candidate to be the source of the disorder. In this regards, it should be noted that the optional module 285 of system 200B may be a computerized module configured and operable for carrying out operation 185 for repeating the determination of the relative phases for several relevant frequency bands and time domain, while possibly receiving as input, or as predetermined data, the relevant bands and time windows, in which to operate.

In some implementations of the present invention the method 100C is implemented for determining a source region sourcing an epileptic seizure in a nervous system. This is based on the assumptions that epileptic seizures start at specific location(s) and spread like waves (i.e., existence of a global temporal 'first') and (2) there are no global temporal 'last' (i.e., specific region(s) in which seizures stop). This is implemented by conducting the operations of method 100C for characterizing the epileptic disorder in the nervous system by modeling n*K high order (e.g. 0>2) temporal hierarchy in networks of m coupled time-signals (e.g. utilizing the method 100A and the operation 140 described above).

In specific tests of the techniques conducted for identifying the source of the epileptic disorder n EEG time-signals obtained from n respective regions/electrodes in the brain were processed. n*K*J temporal hierarchy in networks of 2^nd order (0=2) of m of the time-signals, were modeled utilizing the method 100A and the operation 140 described above. Namely per each candidate signal of the n signals, 2^nd order networks were modeled for k subsets including the candidate signal and this was conducted for J different frequency-bands - time-window portions of the signals).

In other words the operations of method 100C for determining the source of the epileptic seizure were performed on several frequency-bands and several time- windows. As illustrated in Fig. 4, in each frequency band and time-window the implementation of method 100C included the following steps: (1) Random selection of 3 electrodes (i.e. time-signals//zones) and assign one of them to be the reference electrode. (2) Construct m=4 time-signals FC networks (see method 100A) with the three selected electrodes and with a time-signal of another (4^th) electrode, for all electrodes/signals. For each of these networks, determine the hierarchical order (see operation 140 above) to obtain a phase value for the 4^th electrode. This phase value is the (nonlinear) phase-difference of the (4^th) electrode with the reference electrode and gives the electrode temporal rank (order) in the network. (3) Repeat steps (l)-(2) multiple times where in each, a set of n phase-differences are obtained, one for each of the n available electrodes (signals/regions). (4) Sum these phase-differences across all trials. Phase-differences of electrodes whose temporal/hierarchical location within the networks varies, will be suppressed nearly to zero while phase-differences of electrode(s) that correspond to 'global first' will sum to relatively substantial nonzero/ values. The nervous system regions associated with these electrodes (the 'global first' electrodes) are hypothesized to present the epileptogenic zones. Note that phase sigh is a function of the frequency and the time-window therefore, when performing the above processing on several frequency bands and/or several time windows, the absolute phase values are summed in (4) and used to infer epileptogenic zones. Also note that when utilizing different frequency band portions and/or time window portions of the signals, each of the FC and hierarchical networks determined in (2) is determined based on subsets of m=4 corresponding signal portions of the n signals, (whereby here the term corresponding signal portions pertains portions of the m signals which are associated with the same frequency band and same/simultaneous time window).

Figs. 5 and 6A to 6C are graphical illustrations of test study - conducted using the technique of the present invention for identifying epileptogenic zone of a patient. A 35 minutes subdural and intracerebral depth electrodes recording from epileptic patient was used for initial test of the analyses. The data was recorded using 108 contacts. Sixty-four of them were placed on an 8x8 subdural grid on the left fronto-temporal cortex. Other contacts were along six and eight-electrode subdural strips and depth electrodes that were implanted to the hippocampus, the insula, the anterior temporal region, and the subtemporal area. Data were sampled at 2 kHz and was down sample to 256Hz. Some of the electrodes showed noisy signals and were not used in the analysis below. A first aimed of the test was to identify times of seizures and their corresponding frequencies. Based on the assumption that seizures are characterized by signals with high amplitude, all time-signals were decomposed into time-frequency space using wavelet analysis (as described in method 100A above), and their average amplitudes (in each frequency-band) were calculated.

Fig. 5 shows ten graphs of the average wavelet power along time of 10 preselected frequency bands, of iEEG recordings/signals obtained along times extending from before and until after occurrence of epileptic seizures. The recordings/signals where obtained from n electrodes (e.g. 100), and the wavelet power per frequency band in the graphs are averaged across all the electrodes. The X axis in the figure are time, the Y axis is the average wavelet power, and the different graphs pertain to the different frequency bands whereby the average frequencies of the bands are listed on the right and are given in Hz. Each time point of the iEEG signals was the averages across half a second.. From the figure it can be deduced that during recording the patient had 3 seizures, each lasted about one minute. It also suggests that seizure activity appears earlier in time and is of higher power at higher frequencies. This figure was used to obtain the approximate seizures' onset times. Analysis was performed at 15 time-windows of 1 second duration each before these times at each frequency band separately. For the computation approach 500 permutations of groups of m-1 signals with customer design software written in DDL.

From these graphs it appears that epileptic seizures are significant at higher frequencies and therefore in the technique of the present invention may be implemented at even higher frequencies than those used in the test.

Figs. 6A and 6B show the results obtained by the technique of the present invention (method 100C) for determining the source region(s) epileptic seizures of the epileptic seizures. Fig. 6B shows the phase values (averaged across all time-windows and the three seizures) for each of the 10 frequency-bands with frequencies in Hz shown on the left. This '2D' presentation helps to identify electrodes that could infer epileptogenic zone. Fig. 6A shows the absolute value summation of these phases across all frequencies. Absolute values were used since only phase amplitude are count for the analysis. Electrodes number 64, 61 and 53 (all within the 8x8 grid) were significantly above others (phase>3SD). Note however that in Fig. 6B also other electrodes show relatively high phase amplitudes. This is due to the fact that high phase values at only few frequencies might indicate epileptogenic zone in that frequency.

Fig. 6C shows a 2D mage of the open brain of the patient with the electrode grid and location of the depth electrodes. In this image, electrodes marked in yellow correspond to the clinical prediction of epileptogenic zone that was obtained by experience epileptologists using visual inspection of the recordings. Areas marked in white correspond to areas that can't be resected due to their eloquent function. The source regions of the epileptic seizure predicted by the technique of the present invention (method 100C) are marked in blue. As seen our prediction is overlaps with clinical prediction.

Computer Simulations

region-FC analysis. To minimize transient effects, averaging in wavelet space was performed on points 2000-4000 (20 to 40 seconds) and all wavelet frequencies were averaged together.

Fig. 8 illustrates the dependency on Gaussian noise of the 4-region-FC and pairwise coherence phases. The phases were calculated with coupled time-series functions obtained by computer simulations of the Kuramoto model as discussed above. The model was run 100 times for each noise and coupling strength level and the SEs of the phases calculated by 4-region-FC and by pairwise coherences were obtained. The means of these SEs along coupling strength demonstrate lower noise in the 4-region-FC calculations. To this end, Comparison of phase fluctuations between 4-region-FC and pairwise calculations, for coupling constants between 0.5 to 8, shows that the 4-region- FC analysis is less sensitive to Gaussian noise for the high noise levels .

Further to the above described computer simulations, the validity of the technique of the present invention was also demonstrated on resting-state functional MRI data of 40 young healthy subjects to characterize the ventral visual system, motor system and default mode network (DMN). These systems were tested to determine whether they are composed of linear or combined networks, and whether their temporal hierarchy is well-defined. The ventral visual system was selected as a test system since its hierarchy is well known, the motor system was selected to test how recurrent connectivity (i.e., the coexistence of bottom-up and top-down connections) affects the hierarchical architecture, and the DMN was selected as its hierarchy is not known. The ventral visual system is used here as a test system to examine the temporal hierarchy obtained by the new 4-region-FC analysis, since its hierarchy is relatively well established. The organization of the primate ventral visual system is known to have a sequential order of the visual cortical areas from the visual sensory periphery (the primary visual area or calcarine sulcus) to "higher-level" areas involved in abstract aspects of vision, was studied. Using Equation 7, 4-region-FC networks were calculated by means of 3 -seed-Statistical Parametric Maps (hereafter "3-seed-SPM"). 3-seed-SPMs are seed-like SPMs but instead of using a single seed, three seeds are used and 4-region- FC networks are calculated to construct the 3-seed-SPM. For the 4-region-FCs calculations, we used the three time-series functions of the three seeds and a time-series function of a voxel, for all voxels in the brain. In this way, N 4-region-FC networks were used to construct a 3-seed-SPM with N the number of voxels.

presented for wavelet scale 1, with similar results obtained for scale 2. Using Equation 7, the following phases were obtained (mean ± SE):

correspond to time latencies of 0.19 seconds from BA17 to BA18, 0.83 seconds from BA18 to BA37 and 1.33 seconds from BA37 to BA21. These results fit the known temporal organization of the ventral visual system.

Second, 3-seed-SPMs, was calculated to identify functional networks of the ventral visual system in the entire brain. The following regions(seeds) were used: left calcarine sulcus, left cuneus and left fusiform gyrus. Fig. 9 shows the amplitude 3- seeds-SPM and joint-seed-SPM (that is obtained by linear summation of pairwise interactions) for wavelet scale 1 with similar results shown for scale 2. Section (A) of the figure shows the Amplitude 3-seed-SPM, and section (B) of the figure shows the joint-seed-SPM. The amplitude 3-seed-SPM and the joint-seed-SPM are similar to each other with larger connectivity volumes found in the 3-seed-SPM particularly in the occipital and cerebellar regions. For scale 1, the 3-seed-SPM resulted with 26482 significant voxels and the joint-seed-SPM with 18930 significant voxels, 98% of them common to both analyses. These numbers and the visual inspection of the SPMs suggest that the majority of the quadratic networks of the ventral visual system are linear. Note that both analyses identified the sensorimotor cortex to be functionally connected to the ventral visual system.

time-lags of the BOLD signals correspond to the latencies of information flow, the directed hierarchy was obtained for all significant clusters in the brain. As seen, the ventral visual system exhibits a temporal hierarchy generally in line with its known organization.

The motor system is used here as an example for a system that is far from being linear. Functional networks of the motor system were identified by using the following seeds: the left precentral gyrus, left supplementary motor area (SMA) and left thalamus. Fig. 11 presents the amplitude 3-seed-SPM and the joint-seed-SPM for wavelet scale 1. Section (A) shows the amplitude 3-seed-SPM, and section (B) shows joint-seed-SPM. In all scales, the 3-seed-SPMs and the joint-seed-SPMs were markedly different: in scale 1, the 3-seed-SPM resulted with 7080 significant voxels while the joint-seed-SPM resulted with 321 significant voxels, 31% of them common to both analyses. The 3- seed-SPMs included large portions of the precentral gyrus, postcentral gyrus, SMA, thalamus, posterior putamen (in scale 2), and dorsal/posterior insula (scale 2) in both hemispheres. The joint-seed-SPM included smaller parts of the precentral gyrus, postcentral gyrus, thalamus and insula (scale 1) and pre- and post-central regions in both hemispheres in scale 2. The low number of voxels that were simultaneously significant in both; the 3-seed-SPM and the joint-seed-SPM suggests that only few networks in the motor system are linear. Consequently, we suggest that the motor system, as defined by our seeds, is mainly composed of combined quadratic networks.

Fig. 12 shows a phase 3-seed-SPM of the motor system for wavelet scale 1

(0.02-0.03Hz). Section (A) show the phase 3-seed-SPMs for wavelet scales 1 and 2 with the phase of the thalamus as a reference. Section (B) demonstrates the sensitivity of the analysis: by zooming into the thalamus and presenting the phases with finer color bar, we show that the analysis discriminates between different thalamic nuclei each with a different phase, as expected by the anatomy (Wu, Yu et al. 2014). Since both bottom-up and top-down connections are expected in the motor system, a phase relation for all 4- region-FC networks in the motor system cannot be assumed and the system's directed hierarchy cannot be obtained. However, the undirected hierarchy can easily be observed and is of center<→out (i.e., sub-cortex <→cortex) organization: thalamus- SMA- primary motor cortex- sensory motor cortex.