Patent ID: 12248628

DETAILED PRESENTATION OF PARTICULAR EMBODIMENTS

In the following, we will consider a direct neural interface with continuous decoding using a Markov mixture of experts as described in the introductory part.

The electrophysiological signals output from the different electrodes are sampled and assembled by data blocks, each block corresponding to a sliding observation window with width ΔT. For example, the electrophysiological signals are sampled at a frequency of about 590 kHz, and the size of data blocks is 590 samples (ΔT=1 s), and the offset δT from one observation window to the next is 59 samples (δT=100 ms). Each observation window is defined by an observation time (epoch) at which the window in question starts.

The electrophysiological signals are then preprocessed. In particular, this preprocessing may comprise elimination of the average taken on all electrodes, and a time-frequency analysis is then made on each observation window.

The time-frequency analysis can use a decomposition into wavelets, particularly Morlet wavelets.

A frequency smoothing or decimation can also be made on these results of the time-frequency analysis.

Thus, each observation window, or observation instant t, is associated with an order 3 tensor of observation data, resulting in the generation of an order 4 input tensor: the first mode corresponds to successive observation windows, the second mode corresponds to space, in other words the sensors, the third mode corresponds to time within an observation window, in other words the positions of the wavelets, and the fourth mode corresponds to the frequency, in other words to the number of frequency bands used for the decomposition into wavelets on an observation window.

More generally, the input tensor (or the observation tensor) will be order n+1, the first mode always being the mode relative to observation times (epochs). The input tensor (or observation tensor) is denoted X and its dimension is N×I1× . . . ×In.

Similarly, the trajectory of the imagined movement, observed or performed, is described by an output tensor (or command tensor), order m+1, denotedY, with dimension N×J1× . . . ×Jm, of which the first mode corresponds to successive times at which the commands will be applied (in general, this first mode also corresponds to the observation windows), the other modes corresponding to commands on the different effectors or different degrees of freedom of a multi-axis robot.

More precisely, the output tensor supplies N consecutive command data blocks, each of the blocks being used to generate command signals for the different effectors or degrees of freedom. Thus, it will be understood that the dimension of each data block can depend on the envisaged usage case and particularly on the number of degrees of freedom of the effector.

In the following,Xtwill be used to denote the observation tensor at time t. Consequently, this tensor is order n and its dimension is I1× . . . ×In. It takes its values in a spaceX⊂I1× . . . ×Inin whichis the set of reals. Similarly,Ytwill be used to denote the command tensor at time t. Consequently, this output tensor is order m and dimension J1× . . . ×Jm. It takes its values in a spaceY⊂J1× . . . ×Jm.

It will be assumed that the spaceXis formed from the union of a plurality K of regions, not necessarily discontiguous. In other words,

X¯=⋃k=1KXk
in which Xk, k=1, . . . , K are the regions in question.

The direct neural interface is based on a mixture of experts (ME), each expert Ekoperating on the elementary region Xkof the input characteristics space and being capable of predicting the command tensorYtat time t starting from the observation tensor,Xt, when this tensor belongs to Xk. In other words, each region Xkhas a corresponding command tensor prediction modelYtstarting from the observation tensorXt. An expert Ekcan thus be considered as a multi-linear application of XkinY.

It will also be assumed that each expert Ekis associated with a hidden state k of a first order Markov Model (HMM). The different experts are combined by means of combination coefficients dependent on the hidden state at the time of the prediction. Definitively, starting from an input (or observation) tensorXt, the neural interface estimates the output (or command) tensorYt, using:

Y^_t=∑k=1Kγkt(β_k⁢X_t+δ_k)(2)
in whichβk,δkare a prediction coefficients tensor and an expert prediction bias tensor respectively, Ekand δktis the weighting coefficient (also called the gating coefficient) for this expert and at the time t considered. The weighting coefficient γktis simply the conditional probability that the HMM model is in the state k knowing the previous input tensorsX1:t=X1,X2, . . . ,Xt.

The set of coefficients and prediction biases, collectively called prediction parameters of the different experts, is designated by θe={(βk,δk)|k=1, . . . , K}. The set of parameters for the combination of the different experts (including the parameters of the subjacent HMM model) is designated by θg={A,{dk|k=1, . . . , K}, π} in which A is the transition matrix between states, with size K×K, in other words the matrix for which the elements aij(independent of time by assumption of the HMM model) are defined by aij=p(zt=j|zt-1=i) in which zt, represents the state of the model at time t and zt-1represents the state of the model at the preceding time; {dk|k=1, . . . , K} the parameters used to determine the conditional probability of observing the input tensorXtknowing the state zt=k, and π is a vector with size K giving probabilities of occupancy of the different states at the initial time, in other words πi=p(zt=i)t=0.

FIG.2schematically represents a direct neural interface using a Markov mixture of experts according to one embodiment of the invention.

The electrophysiological signals output from the different electrodes205are processed in the processing module210. This processing module performs sampling, optionally eliminates the average taken on all electrodes, and a time-frequency analysis is then made on each observation window. The time-frequency analysis can possibly be followed by frequency smoothing or decimation, as indicated above. Consequently, the processing module210outputs an input tensorXtat each observation window, characterised by an observation time t.

The direct neural interface also makes use of a hidden state machine based on an HMM model,240, that can be in K possible states, and a plurality of experts,220, each expert being associated with a hidden state.

Each expert Ekmakes an estimateYtkof the output vectorYtand supplies it to a combination module,230, that mixes the estimates according to expression (2). The combination coefficients (gating coefficients) are determined in the estimating module,250, by γkt=p(zt=k|X1:t), k=1, . . . , K, in which ztrepresents the state of the model at time t.

The estimateŶtgives the different kinematic movement parameters such as the position and/or the speed and/or the acceleration along the required trajectory.

The neural interface is calibrated firstly during an initialisation phase and then, during subsequent calibration phases, at given time intervals (on-line calibration).

These subsequent calibration phases are necessary due to the lack of stationarity of signals generated in the human brain. They can be done at regular intervals along the trajectory or successive trajectories to be executed.

A calibration phase u uses a plurality of input data blocks corresponding to a plurality ΔL of successive times output from the observation times sequence.Xuwill denote the tensor with dimension ΔL×I1× . . . ×Incontaining this plurality of input blocks. It also makes use of a plurality of output blocks giving kinematic set parameters, taken at the same times.Yuwill denote the set tensor with dimension ΔL×J1× . . . ×Jmcontaining this plurality of input blocks.

The switches261,262are in the closed position during a calibration phase. The calibration module270receives the input tensorXu, the set tensorYu, and the matrix Zuwith size K×ΔL, of which each of the columns represents the state of the machine at the time considered. More precisely, if at timethe machine is in state z=k, all elements in thethcolumn of the matrix Zuare null except for the element on row k that is equal to 1.

The calibration module usesXu,Yu,Zuto calculate the updated parameter sets θeand θg.

At the end of this calibration phase, the switches261,262are open and the updated parameter sets θeand θgare output to a plurality of experts,220, and to the machine240, respectively.

Innovatively, the predictive model of each expert Ekis trained using an REW-NPLS (Recursive Exponentially Weighted N-way Partial Least Squares) regression with a forget factor λk(0<λk<1) on calibration data blocks,Xu,Yu, corresponding to the hidden state zu=k, denotedXukandYuk. More precisely, for each expert Ekand for each calibration phase u, the data blocksXukandYukare determined so that it can be trained in a supervised manner.

This is done as described in application FR-A-3061318, namely:

In a first step, the tensorsXukandYukare normalised starting from a number of observations said to be effective Nuk=λkNu-1k+Nukin which Nukis the number of data blocks corresponding to the hidden state zu=k during the calibration phase u. The weighted sum s(Xuk) and the weighted quadratic sum sq(Xuk) are calculated:

s⁡(X_uk)=λk⁢∑ℓ=1Nu-1k⁢xu-1;ℓ,i1,⁢…⁢⁢ink+∑ℓ=1Nuk⁢xu;ℓ,i1,⁢…⁢⁢ink(3⁢-⁢1)sq⁡(X_uk)=λk⁢∑ℓ=1Nu-1⁢(xu-1;ℓ,i1,⁢…⁢⁢ink)2+∑ℓ=1Nuk⁢(xu;ℓ,i1,⁢…⁢⁢ink)2(3⁢-⁢2)in whichandrepresent the elements of tensorsXu-1kandXuk, respectively. The average value

μ⁡(X_uk)=s⁡(X_uk)Nuk
and the standard deviation

σ⁡(X¯uk)=s⁢q⁡(X¯uk)-(μ⁡(X¯uk))2Nuk
are then deduced. The elements of the centred and normalised input tensor,{tilde over (X)}uk, are:

x~u;ℓ,i1,⁢…⁢⁢ink=xu;ℓ,i1,⁢…⁢⁢ink-μ⁡(X_uk)σ⁡(X_uk)(4)

Similarly, the weighted sum s(Yuk) and the weighted quadratic sum sq(Yuk) are calculated:

s⁡(Y¯uk)=λk⁢∑ℓ=1Nu-1k⁢yu-1;ℓ,j1,⁢…⁢⁢jmk+∑ℓ=1Nuk⁢yu;ℓ,j1,⁢…⁢⁢jmk(5⁢-⁢1)sq⁡(Y¯uk)=λk⁢∑ℓ=1Nu-1k⁢(yu-1;ℓ,j1,⁢…⁢⁢jmk)2+∑ℓ=1Nuk⁢(yu;ℓ,j1,⁢…⁢⁢jmk)2(5⁢-⁢2)

The average value

μ⁡(Y_uk)=s⁡(Y¯uk)Nuk
and the standard deviation

σ⁡(Y¯uk)=sq⁡(Y¯uk)-(μ⁡(Y¯uk))2Nuk
are then deduced, followed by the elements of the centred and normalised output tensor,{tilde over (Y)}uk:

y~u;ℓ,j1,⁢…⁢⁢jmk=yu;ℓ,j1,⁢…⁢⁢jmk-μ⁡(Y¯uk)σ⁡(Y¯uk)(6)

The covariance tensor and the cross-covariance tensor are defined during the calibration phase u, starting from the centred and normalised tensors{tilde over (X)}ukand{tilde over (Y)}uk:

cov⁡(X~_uk,X~_uk)=X~_uk×1X~_uk(7-1)cov⁢(X~_uk,Y~_uk)=X~_uk×1Y~_uk(7-2)
in which x1designates the tensor product according to the first mode (temporal mode with index) as described in the above-mentioned application FR-A-3061318.

These covariance tensors are modified taking account of the covariance tensors of the previous calibration phase, by weighting them with the forget factor λk:
cov({tilde over (X)}uk,{tilde over (X)}uk)=λkcov({tilde over (X)}u-1k,{tilde over (X)}u-1k)+{tilde over (X)}uk×1{tilde over (X)}uk(8-1)
cov({tilde over (X)}uk,{tilde over (Y)}uk)=λkcov({tilde over (X)}u-1k,{tilde over (Y)}u-1k)+{tilde over (X)}uk×1{tilde over (Y)}uk(8-2)

This calculation takes account of the data from a previous calibration phase, with a forget factor, to update the covariance tensors. The forget factors relative to the different experts can be chosen to be identical λk=λ, k=1, . . . , K, in which λ is then the common forget factor. Alternatively, they can be chosen to be distinct so as to offer more prediction flexibility to the different experts.

Starting from the covariance tensors cov({tilde over (X)}uk,{tilde over (X)}uk) and cov({tilde over (X)}uk,{tilde over (Y)}uk) working iteratively on rank fk=1, . . . , F in the latent variables space (in which F is a maximum predetermined number of dimensions in the latent variables space), we obtain a set of projection vectors {wu,1k,fk, . . . , wu,nk,fk}fk=1Fwith dimensions I1, . . . , Inrespectively, a set of prediction coefficient tensors, {βuk,fk}fk=1F, and a prediction bias sensor, {δuk,fk}fk=1F, using the REW-NPLS method described in the above-mentioned application FR-A-3061318.

Starting from the covariance tensors cov({tilde over (X)}uk,{tilde over (X)}uk) and cov({tilde over (X)}uk,{tilde over (Y)}uk) working iteratively on rank fk=1, . . . , F in the latent variables space (in which F is a maximum predetermined number of dimensions in the latent variables space), we obtain a set of projection vectors {wu,1k,fk, . . . , wu,nk,fk}fk=1Fwith dimensions I1, . . . , Inrespectively, a set of prediction coefficient tensors, {βuk,fk}fk=1F, and a prediction bias tensor, {δuk,fk}fk=1F, using the REW-NPLS method described in the above-mentioned application FR-A-3061318.

The optimum rank, foptk, for each expert Ek, during a calibration phase u is determined by calculating the prediction errorsYukstarting fromXukusing the prediction coefficient tensors, {βuk,fk}fk=1F, and the prediction bias tensors, {δuk,fk}fk=1Fobtained during the previous calibration phase:
erruk,fk=λkerru-1k,fk+∥Yuk−Ŷu-1k,fk∥  (9)
in whichŶk,fkis the estimate ofYukobtained from the predictionβu-1k,fkand prediction bias coefficientsδu-1k,fkfrom the previous calibration phase.

The rank foptkleading to the minimum prediction error is then selected:

fo⁢p⁢tk=arg⁢⁢minfk=1,⁢…⁢,F⁢(erruk,fk)(10)
and for the calibration phase, the prediction coefficients tensor and the prediction bias tensor corresponding to this optimal rank are selected, in other words the update is made at the end of the calibration phase u:
βk=βuk,foptk(11-1)
δk=δuk,foptk(11-2)
for each of the experts Ek, k=1, . . . , K, independently. The result obtained is thus an estimate of the set of prediction parameters θemaking use of the REW-NPLS regression method applied to each of the experts.

Innovatively, the combination coefficients γkt, k=1, . . . , K are also estimated using an REW-NPLS regression method adapted to discrete decoding, as described later.

During a calibration phase, the elements aijof the transition matrix A are firstly updated as follows:

auij=λ⁢au-1ij+vuijK(12)
in which auij(or au-1ij) are elements of the transition matrix, estimated during the calibration phase u (or u−1) and is the number of transitions from state i to state j during the calibration phase u. This number of transitions is obtained starting from the matrix Zu. The lines in the transition matrix are then normalised such the sum of the probabilities of the transition from an initial state is equal to 1.

The transition matrix A thus updated is supplied by the transition matrix calculation module,260, to the state estimating module,240, of the HMM model.

zldefines a vector with size K for which the elements represent whether or not the machine is in state k=1, . . . , K at time t. Thus, if the machine is in state k at time t, only the kth element of ztwill be equal to 1 and the other elements will be null. ztcan be considered like a process with discrete values that vary in time according to a multilinear predictive model, the input variable of which is the observation tensorXt. This predictive model is trained under supervision during each calibration phase u, using the same principle as the predictive model of each expert. In other words:
{circumflex over (z)}t=BXt+b(13)
in which {circumflex over (z)}texpresses the probabilities {circumflex over (z)}k,tthat the machine is in state k at time t,Bis a prediction coefficients tensor with size K×I1× . . . ×Inand b is a bias vector with size K.

The tensorBand the vector b are updated during each calibration phase u by means of an REW-NPLS regression with a forget factor i λ, (0<λ<1) based on observationsXuand Zu(matrix with size K×ΔL). This forget factor may be identical to the common forget factor, when such a common forget factor is used to estimate prediction parameters for the different experts.

Working iteratively on rank f=1, . . . , F of the latent variables space starting from the covariance tensor cov({tilde over (X)}u,{tilde over (X)}u) in which{tilde over (X)}uis a centred and normalised version ofxu, and the cross-covariance tensor cov({tilde over (X)}u, Zu). These covariance and cross-covariance tensors are modified by:
cov({tilde over (X)}u,{tilde over (X)}u)=λ·cov({tilde over (X)}u,{tilde over (X)}u)+{tilde over (X)}u×1{tilde over (X)}u(14-1)
cov({tilde over (X)}u,Zu)=λ·cov({tilde over (X)}u,Zu)+{tilde over (X)}u×1Zu(14-2)

The REW-NPLS method provides a set of projection vectors {wu,1, . . . , wu,n}f-1F, with dimensions I1, . . . , Inrespectively, a set of prediction coefficient tensors, {βuf}f-1F, and a set of prediction bias vectors, {buf}f-1F. Rank foptleading to the minimum prediction error is selected, then the tensorBand the vector b are updated by:
B=Bufopt(15-1)
b=bufopt(15-2)

The tensorBand the vector b are supplied to the mixture coefficients estimating module,250.

Starting from elements {circumflex over (z)}k,tof the predicted vector {circumflex over (z)}t, the state estimating module can calculate the conditional probabilities using the softmax function:

p⁡(zt=k|X¯t)=exp⁡(z^k,t)∑i=1K⁢exp⁡(z^i,t)(16)

This expression gives the probability that the machine is in state k at time t, knowing the input data at this time.

The mixing coefficients estimating module,250, obtains these coefficients from Bayes rule:

γkt=p⁡(zt=k|X_1:t,)=p⁡(zt=k,X_1:t)p⁡(X_1:t)=p⁡(zt=k,X_1:t)∑i=1K⁢p⁡(zt=i,X_1:t)(17)
in which p(zt=i,X1:t) are obtained using the forward algorithm, namely:

p⁡(zt=i,X_1:t)=p⁡(X_t|zt=i)·∑j=1K⁢aij⁢γkt-1(18)

The mixture coefficients can be obtained by recurrence by combining expressions (17) and (18):

γkt=p⁡(X_t|zt=k)·∑j=1K⁢akj⁢γkt-1∑i=1K⁢(p⁡(X_t|zt=i)·∑j=1K⁢aij⁢γit-1)(19)

The conditional emission probabilities, p(Xt|zt=i) can be obtained using Bayes rule starting from a posteriori conditional probabilities p(zt=i|Xt) and a priori probabilities p(zl=i|Xt):

p⁡(X_t|zt=i)=p⁡(zt=i|X_t)⁢p⁡(X_t)p⁡(zt=i)∝p⁡(zt=i|X_t)p⁡(zt=i)(20)
the term p(Xt) being a multiplication coefficient common to all states.

The probabilities of occupancy of the different states p(zt=i) at time t are calculated from the initial probabilities of occupancy π and the transition matrix A.

FIG.3is a flowchart representing a method of iteratively calibrating a direct neural interface using a Markov mixture of experts according to one embodiment of the invention;

Said calibration method works during calibration phases planned at determined instants, for example at regular time intervals throughout the trajectory. Each calibration phase u uses a plurality of input data blocks corresponding to a plurality ΔL of successive times in the observation times sequence. It also makes use of a plurality of output blocks giving kinematic set parameters at the same times.

The plurality of input data blocks for the calibration phase u is represented by the input tensorXuand the plurality of output blocks during this same calibration phase is represented by the output tensorYu.

It is also assumed that the states of the HMM machine are known during the calibration phase. The different states of the HMM machine can relate to different elements to be controlled to make the trajectory, for example different members of an exoskeleton. Furthermore, for each of these elements, different states can be considered for example such as an active state when the patient controls this element, an idle state when he does not control it, and possibly a preparation state immediately preceding the active state and immediately following the idle state, during which the characteristics of the neural signals are not the same as in the idle state nor the active state for this element. Thus, when the direct neural interface must control P elements to make the trajectory, the HMM machine will comprise 2Por even 3Pstates (if preparatory states are envisaged). The states of the HMM machine during the calibration phase can be represented by a binary matrix Zuwith size K×ΔL, all values in each column in Zubeing 0 except for one element equal to 1 indicating the state of the machine at the corresponding time.

The input tensor,Xu, the output tensor,Yu, the states matrix, Zu, are given in310for the calibration phase u.

In step320, the matrix Zuis used as a starting point to determine the tensorsXukandYukformed by the input data blocks and the output data blocks related to the state k=1, . . . , K, respectively. These tensors are extracted from the input tensor,Xu, and the output tensor,Yu.

In step330, for each expert Ei, k=1, . . . , K the centred and normalised tensors{tilde over (X)}ukand{tilde over (Y)}ukare calculated and are used to deduce the covariance tensor{tilde over (X)}uk, and the cross-covariance tensor of{tilde over (X)}ukand{tilde over (Y)}uk. These tensors are modified by adding the covariance tensor of{tilde over (X)}u-1kand the cross-covariance tensor respectively, obtained in the previous calibration phase, weighted by a forget factor λk. These tensors thus modified will be used as covariance and cross-covariance tensors during the next calibration phase.

In step340, the prediction parameter tensorsβk,δkof each expert Ekare updated using an REW-NPLS regression, starting from the covariance and cross-covariance tensors modified in the previous step.

After step340, we have an updated version of the set θeof expert prediction parameters K.

In step350, the elements of the transition matrix A are advantageously updated starting from the number of transitions between successive states observed during the calibration phase u. The number of transitions between successive states is obtained starting from the matrix Zu. The matrix A is then modified by adding to it the matrix A obtained during the previous iteration, weighted by a forget factor λ.

In step360, the centred and normalised tensor{tilde over (X)}uis added, followed by the covariance tensor of{tilde over (X)}uand the cross-covariance tensor of{tilde over (X)}uand Zu. These tensors are modified by adding the covariance tensor of{tilde over (X)}u-1kand the cross-covariance tensor{tilde over (X)}u-1kand Zu-1respectively, obtained in the previous calibration, weighted by the forget factor λ. These tensors thus modified will be used as covariance and cross-covariance tensors during the next calibration phase.

In step370, a multi-linear predictive model is trained giving a state vector of the HMM machine as a function of the input tensor, {circumflex over (z)}t=BXt+b, the components of the state vector {circumflex over (z)}tproviding probabilities that the machine is in the different states k=1, . . . , K respectively at time t. More precisely, the prediction parameter tensorB, and the bias vector b are updated using an REW-NPLS regression, starting from the covariance and cross-covariance tensors modified in the previous step.

The components of {circumflex over (z)}tcalculated during an operational phase are then used at each observation time t to obtain mixture coefficients γkt=p(zt=k|X1:t) of the different experts.