Computer-implemented method for inferring a conduction velocity of a subject's brain

Abstract

A way to infer a conduction velocity of a subject's brain, from a recorded description of the frequency content of the activities of the subject's brain, comprises the steps of: (i) customizing a digital brain model with neuronal activity data recorded on the subject, the digital brain model being capable of simulating realistic personalized activities of the subject's brain in the form of time-series, (ii) inputting n plausible values of average conduction velocity into the customized digital brain model to obtain time-series in order to characterize n descriptions of the frequency content of the activities of the simulated brain, (iii) establishing a backward function from the calculated descriptions to the inputted values, (iv) applying the backward function to the recorded description of the frequency content of the activities of the subject's brain to obtain the average conduction velocity of the subject's brain.

Description

DescriptionTitle: Computer-implemented method for inferring a conduction velocity of a subject's brainTechnical field

[0001] The technical field of the invention is brain function simulation, and more specifically the assessment of conduction velocity and its applications to myelin-related diseases, such as multiple sclerosis.State of the art[2] WO 2014 / 020428 discloses a method and apparatus for simulating brain function. It makes available a brain model that is configured thanks to measured neuroimaging data. The content of this document is incorporated by reference in the present description.[3] US-2023 / 0022257 discloses a data processing system configured to compute, based at least in part on the medical-imaging data, a set of white-matter disease biomarkers for different neurological anatomical regions. This system is used in relation to the description and prediction of neuro-degenerative disease.[4] The system is able to provide feedback information such as diagnostic information, information associated with disease progression, information regarding efficacy of a treatment, or a treatment recommendation.Technical problem to be solved[5] In order to improve a diagnostic, the skilled person is always looking for more accuracy. However, increase of accuracy usually goes with increase of resources and computing power. One ultimate purpose of the invention is to provide an alternative way to predict the outcome of a neuro-degenerative disease or at least to monitor the progression of such a disease. Along the way, the invention provides other objects.Disclosure of the invention[6] An object of the invention is a computer-implemented method for inferring an average conduction velocity of a subject's brain from a recorded description of the frequency content of the activities of the subject's brain, characterized in that it comprises the steps of:(i) customizing a digital brain model with neuronal activity data recorded on the subject, the digital brain model being capable of simulating realistic personalized activities of the subject's brain in the form of time-series;(ii) inputting n plausible values of average conduction velocity into the customized digital brain model to obtain time-series in order to characterize n descriptions of the frequency content of the activities of the simulated brain;(iii) establishing a backward function from the calculated descriptions to the inputted values;(iv) applying the backward function to the recorded description of the frequency content of the activities of the subject's brain to obtain the average conduction velocity of the subject's brain.[7] One advantage of this method is to provide a mechanistic link that makes it possible to access the value of conduction velocity of a subject's brain. This indicator is a good candidate to measure sub-clinical damage to the neural structures, in particular the myelin sheaths.[8] In the present description, the following terms are used with the following meanings:[9] -subject: mammal, preferably human being;

[0010] -neuronal activity data: data that reflects ongoing activity from within neural tissues;

[0011] -brain regions: functionally or structurally homogeneous brain areas that are conceptualized as separated entities;

[0012] -frequency content of the activities of a brain: frequencies contained in timeseries that have been generated by neuronal activities;

[0013] -plausible values: values that are considered plausible according to the previous knowledge of the field;

[0014] -invertible function: a function f such that another function f-1 exists as to undo the operation of f;

[0015] -forward function: a map from the velocities to the descriptions;

[0016] -backward function: a map from the descriptions to the velocities;

[0017] -global coupling value: a number that quantifies the importance of the long-range connections in the brain model.

[0018] According to a particular embodiment of the method, the digital brain model is a network of a plurality of brain regions, each brain region being modelized by a local area model, the local area models being connected to each other in the network configured according to a fractography carried out on the subject.

[0019] One advantage of the network of brain regions is the possibly light implementation of the method, in comparison to resource-consuming models that can only increase the definition by scaling up the size of datasets.

[0020] According to a particular embodiment of the method, any inputted average conduction velocity value is combined with any global coupling value.

[0021] One advantage of these parameters is their availability according to state-of-art standards.

[0022] According to a particular embodiment of the method, the digital brain model is based on a 3D-map of the cerebral cortex.

[0023] One advantage of the 3D-map model is its lower dependency on the definition of the brain regions.

[0024] According to a particular embodiment of the method, the backward function is established using an artificial neural network.

[0025] One advantage here is the ability of the approach to set up a relationship between the input values and the expected output values, even when a closed-form solution cannot render the service.

[0026] According to a particular embodiment of the method, the backward function is obtained by first establishing an invertible forward function from the inputted values to the calculated descriptions, then inverting the forward function from the descriptions to the conduction velocities.

[0027] One advantage of going through an invertible function is the constructional approach: since the dataset inputs the descriptions and outputs the conductivity, and possibly coupling parameter, it may become easy to first obtain an invertible function, especially thanks to a neural network, then to invert it to obtain the forward function.

[0028] According to a particular embodiment of the method, the description of the frequency content of the activities of a brain is a power spectra description.

[0029] One advantage of the PSD is that it is well-known and available instrument that can be calculated on a number of types of neural data, and on scanners from any vendor. This facilitates one implementation of the invention.

[0030] According to a particular embodiment of the method, the neuronal activity data contains non-invasive magneto-encephalography data or electro-encephalography data.

[0031] According to a particular embodiment of the method, in step, 15,000 to 25,000 plausible values, preferably 20,000, are inputted into the customized digital model.

[0032] One advantage of such figures is to provide a reasonable definition while requiring a reasonable computing power. As explained in the example, variations of these figures increase the definition, but not drastically. These figures are a good compromise.

[0033] According to a particular embodiment of the method, the plausible values are randomly set in a range of 0.5ms to 20ms.

[0034] One advantage of this rather wide range is to cover all the known possible values that may need to be inferred when implementing the method.

[0035] Another object of the invention is a computer program product, characterized in that it comprises instructions which, when the program is executed by a computer, cause the computer to perform the method.

[0036] Another object of the invention is a data processing system, characterized in that it comprises a processor configured to perform the method.

[0037] Another object of the invention is a computer-readable storage medium, characterized in that it comprises instructions which, when executed by a computer, cause the computer to perform the method.

[0038] Another object of the invention is a method for obtaining an indicator of the progression of damage to the myelin sheath of a subject's neurons, comprising the steps of:- comparing an earlier inferred conduction velocity of the subject's brain to a conduction velocity of the subject's brain at a later time point,- providing an indicator of progression as an increasing function of the comparison between the consecutive timepoints, characterized in that the earlier and later inferred conduction velocities are obtained by carrying out the method for inferring a conduction velocity. The comparison may be implemented as a difference.

[0039] According to a specific embodiment, the indicator of progression is set such as an earlier inferred conduction velocity higher than a later inferred conduction velocity is indicative of a poor indicator of the progression of damage to the myelin sheath, whereas an earlier inferred conduction velocity lower or equal than a later inferred conduction velocity is indicative of a good indicator of the progression of damage to the myelin sheath.

[0040] By “a poor indicator of the progression of the damage to the myelin sheath” or “the indicator of progression is high”, it is meant that the myelin sheath is getting damaged over time. Therefore, there is still a progression of the damage to myelin sheath over time.

[0041] By “a good indicator of the progression of the damage to myelin sheath” or “the indicator of progression is low”, it is meant that there is no additional damage to the myelin sheath over time and / or that the myelin sheath is getting resorbed over time. Therefore, there is at least no progression of the damage to myelin sheath over time or even a reduction of the damage of the myelin sheath.

[0042] One advantage of this method is to make use of the conduction velocity as a sensitive maker of the myelin damage, as opposed to solely relying on structural damage measured by structural MRI.

[0043] Another object of the invention is a method for monitoring the efficacy of a medical treatment for multiple sclerosis, comprising the steps of:- obtaining an indicator of the progression of damage to the myelin sheath from a recorded description of the frequency content of the activities of a subject's brain;- indicating that the treatment is efficient when the indicator of progression is low or indicating a lack of efficacy of the therapy when the indicator of progression is high, characterized in that obtaining the indicator results from carrying out the above-described method for obtaining an indicator.

[0044] One advantage of this method is to reverse the usual approach: instead of monitoring structural myelin sheath lesions, a more relevant criteria, conduction velocity, is used to monitor the myelin sheath condition directly from its functioning, in particular in the context of a medical treatment.

[0045] Another object of the invention is an in-silico use of the method for inferring a conduction velocity to quantify and / or assess and / or predict the severity of the disability of a subject with multiple sclerosis, from a recorded description of the frequency content of the activities of said subject's brain.

[0046] Another object of the invention is a method to quantify and / or assess and / or predict the severity of the disability of a subject with multiple sclerosis, using the inferred average conduction velocity of the subject's brain obtained from an in-silico use of the method for inferring a conduction velocity.

[0047] Another object of the invention is a method to detect functional alterations in otherwise normally appearing myelin in the brain of a subject, as demonstrated by a slowing over time of the inferred average conduction velocity of a subject's brain obtained from an in- silico use of the method for inferring a conduction velocity.

[0048] Another object of the invention is a method to detect damage to the myelin sheath of a subject's neurons comprising:-inferring the average conduction velocity of the subject's brain using the inferred average conduction velocity of the subject's brain obtained from an in-silico use of the method for inferring a conduction velocity,-concluding that the subject has a functional damage of the myelin sheath utilizing the inferred conduction velocities; or concluding that the functionality of the myelin sheath of the subject's neurons is unaffected.

[0049] Another object of the invention is a method to monitor the efficacy of a medical treatment for multiple sclerosis from the evolution over time of the inferred average conduction velocity of a subject's brain obtained from an in-silico use of the method for inferring a conduction velocity.

[0050] Another object of the invention is an in-silico use of the method for inferring a conduction velocity to monitor multiple sclerosis by detecting a decline over time of an inferred average conduction velocity of a subject's brain.

[0051] Medical treatments pointed by the invention are typically aimed at preventing the immune system of the patient from attacking its myelin. However, sometimes myelin is beingdamaged although it appears normal at MRI imaging. As such, the rationale of the proposed methodology lies in the fact that spotting subtle alterations of the functionality of the myelin sheath might help assessing the efficacy of the therapy in preventing myelin damage. Examples of such therapies might be:

[0052] • Interferon-beta drugs

[0053] • Glatiramer acetate

[0054] • Dimethyl fumarate

[0055] • Fingolimod

[0056] • Teriflunomide

[0057] • Natalizumab

[0058] • Alemtuzumab

[0059] • Ocrelizumab

[0060] • Cladribine

[0061] • Siponimod

[0062] • OzanimodBrief description of the figures

[0063] The invention will be better understood by reading the description of the following non-limiting examples, illustrated by the following figures:-figure 1 shows an overall pipeline illustrating the steps of the method according to an embodiment of the invention,-figure 2 shows an estimated posterior distribution of the brain parameters velocity V and global coupling parameter G,-figure 3 is a set of graphs evidencing the robustness of the model,-figure 4 is a pipeline illustrating the validation of the model on a cohort of subjects,-figure 5 is a set of graphs detailing the multilinear regression of figure 4,-figure 6 is a bloc diagram showing the steps of the method corresponding to figure 1 , -figure 7 is a bloc diagram showing the steps of a method for obtaining an indicator of the progression of damage to the myelin sheath of a subject's neurons,-figure 8 is a bloc diagram showing the steps of a method for monitoring the efficacy of a medical treatment for multiple sclerosis.Description of at least one way of carrying out the invention

[0064] First, the method 100 for inferring a conduction velocity of a subject's brain, implemented by a computer, will be described.Setting up data

[0065] An MRI recording has been previously performed on the subject with an RMI scanner. Sequences of the recording include echo-planar imaging for DTI reconstruction - TR / TE 12,000 / 95.5 ms, voxel 0.94x0.94x2.5 mm3, 32 equally spaced diffusion-sensitizing directions, 5 B0 volumes - and 3D-FLAIR volume for WM lesion segmentation - TR / TE / TI 7000 / 145 / 1919 ms, echo train length 170, 212 sagittal partitions, voxel size 0.52x0.52x0.80 mm3.

[0066] A processing of the diffusion MRI data is carried out using the software modules provided in the FMRIB software library - FSL, http: / / fsl.fmrib.ox.ac.uk / fsl. The diffusion MRI dataset is corrected for head movements and eddy current distortions using the “eddy_correct” routine, rotating the diffusion sensitizing gradient directions accordingly, and a brain mask is obtained from the B0 images using the Brain Extraction Tool routine.

[0067] A diffusion-tensor model is fitted at each voxel, and fiber tracts are generated over the whole brain by deterministic fractography using the Fiber Assignment by Continuous Tracking - FACT - algorithm implemented in Diffusion Toolkit - angle threshold 45°, spline- filtered, masking by the FA maps thresholded at 0.2. A cortical study-specific Region of Interest - ROI - dataset is obtained by masking the ROIs available in an MNI space-defined volumetric version of the Desikan-Killiany-Tourville - DKT - ROI atlas using the GM tissue probability map available in SPM - thresholded at 0.2. This is done to limit the endpoints of the fibers to cortical and adjacent juxtacortical white matter voxels in the subsequent ROI-based analysis of the fractography data.

[0068] The Fiber Assignment volume is spatially normalized to the Fiber Assignment template provided by FSL using SPM12, and the resulting normalization matrices are inverted and applied to the ROI set.

[0069] A source reconstruction is performed, also using the Fieldtrip Toolbox. In brief, the signal time series are reconstructed using 84 ROIs based on the DKT atlas. The reconstruction takes place using the volume conduction model proposed by Nolte. Based on the native MRIs of each subject, the linearly constrained minimum variance - LCMV - beamformer is applied to reconstruct the signal sources based on the centroids of each ROI.Step i: Customization 110 of a digital brain model with neuronal activity data recorded on the subject.

[0070] Each ROI is conceptualized as an autonomous Stuart-Landau oscillator. Each ROI is regarded as a local, feedback inhibition driven E-l unit, oscillating in the gamma band.

[0071] The Stuart-Landau is set as to show dampened oscillations (a<0).

[0072] Oscillators are connected to each other via white-matter, with coupling strength specified by the subject-specific number of DTI fibers across any tracts.

[0073] Adjacency matrices are scaled by a global coupling parameter (G).

[0074] Coupling between ROIs is subject to finite conduction times, estimated by dividing inter-node Euclidean distances by an average conduction velocity,[Equation 1]The activity of each ROI is given as:[Equation 2]where Z is a complex variable and Re[Z(t)] is the corresponding time-series.

[0075] Accordingly, the natural frequency of each ROI is set to a>j = a>0= 2 * Ti * 40.

[0076] To capture the damped nature of local oscillations, a is set as -5.

[0077] Noise, with amplitude p=1 is added to each oscillator to model stochastic fluctuations.

[0078] This digital model is now ready to generate realistic personalized activities of the subject's brain in the form of time-series.Step ii: Input 120 of plausible values of average conduction velocity into the customized digital brain model

[0079] We use a simulation budget size of 20000 plausible values of global coupling parameter and average conduction velocity. The latter is optional according to the invention, but used in this example. As a result, 20000 descriptions of the frequency content of the activities of the simulated brain are obtained in response to the 20000 plausible values.

[0080] Since the purpose is to model alpha oscillations, only the regime spanned byG = 800 - 1600 and mV = 1 - 20 — s is studied. The parameters are drawn randomly from a uniform prior in the ranges: Ge [800,1600], and Ve [1 ,20],

[0081] Simulations are run for 100s by numerically integrating the system of equations (Equation 2), through the Euler-Murayama method.

[0082] The model simulation and parameter estimation are performed on a Linux machine with Intel Core i9-10900 CPU 2.80GHz and 32 GB of RAM.

[0083] From the realistic personalized activities, one calculates the median PSD, peak alpha frequency, peak alpha spectral density and total alpha power, i.e. the area under thespectral density curve between 8-13 Hz. These three values are also named “amplitude”, “median frequency”, and the “total alpha power”. They are kept as a description of the frequency content of the realistic personalized activities of the simulated brain. When the customized digital brain model is run, these three values are estimated from the personalized synthetic data and constitute a description of the frequency content of the realistic personalized activities of the simulated brain. In other words, by transforming the observation from the timedomain to the frequency-domain, we use the summary statistics of PSD: i.e. amplitude, median frequency, and the total alpha power, as the data feature for the training and the inference steps i.e. as descriptions of the frequency content of the activities of the simulated brain.Step iii: Establishment 130 of a backward function from the calculated descriptions to the inputted values

[0084] First, we build an invertible function that, given any particular average conduction velocity and global coupling parameter, provides us with the expected peak-frequency, peakamplitude, and total power in the alpha band. This is done via a Masked Autoregressive Flow - MAF -, which supports invertible nonlinear transformations. The MAF is trained with the budget of 20000 plausible values of global coupling parameter and average conduction velocity and 20000 descriptions of the frequency content of the activities of the simulated brain, i.e. peak-frequency, peak-amplitude, and total power in the alpha band. The PyTorch-based SBI package is used. More in detail, given data y and model parameters 0, Bayes rule defines the posterior distribution as p(0|y) oc p(y|0)p(0), where p(0) is the prior probability distribution over a plausible range of values, and p(y 10) is the likelihood distribution, which quantifies the probability that the observation y is generated by a specific set of 0. Here, e={G, n, with G defined as the global coupling and V as the average conduction velocities across the white matter bundles, whereas y represents the PSD features - peak frequency, peak amplitude and total power in the alpha band, observed in the MEG data of the subject. The key challenge to efficiently perform Bayesian inference lies in the evaluation of the likelihood function p(y|0), given the samples from the prior distribution. In other words, one needs the calculation of intractable integrals to evaluate the likelihood of observing y given a certain value of 0 (that is, the delays and coupling). Then the likelihood function is taken as the forward function.

[0085] The traditional approaches for SBI suffer from the curse of dimensionality and are sensitive to the ad-hoc choices (i.e., rejection thresholds, distance functions, and summary statistics), which significantly affects both the computational efficiency and the accuracy; on the other hand, the SBI with the deep neural density estimator used in this study transforms asimple distribution conditioned on the data features to obtain the full probability distribution of the target parameter, while dealing with non-linear and high-dimensional latent spaces and highly structured data (without the sensitivity to the tolerance level in the accepted / rejected parameter setting).

[0086] The function thus obtained is smooth and invertible.

[0087] Then, this function is inverted to obtain a backward function from the calculated descriptions to the inputted values, that makes it possible to estimate the most likely individual velocities and coupling parameter, given the observed individual spectral features.Step iv: Application 140 of the forward function to the recorded description of the frequency content of the activities of the subject's brain

[0088] An MEG recording is composed of eyes-closed resting-state segments of 3'30” each. Four anatomical coils are applied on the head of each participant and their position is recorded along with the position of four head anatomical points, to identify the position of the head during the recording.

[0089] MEG processing is performed using the Fieldtrip Toolbox. Eye blinking, if present, and heart activity are recorded through electro-oculogram - EOG - and electrocardiogram - ECG -, to identify physiological artifacts. An expert rater checks for noisy signals and removes them. An anti-alias filter is applied to the MEG signals, acquired at 1024 Hz, before being filtered with a fourth order Butterworth HR band-pass filter of 0.5-48 Hz. Principal component analysis and supervised independent component analysis are used to remove the environmental noise and the physiological artifacts, respectively.

[0090] The empirical PSD is observed on the processed MEG recording and the three values: peak frequency, peak amplitude and total power in the alpha band, are extracted from the Empirical PSD of the subject, as a description of the frequency content of the activities of the subject's brain.

[0091] Then, the continuous function is used to estimate the most likely velocity from said description of the frequency content.

[0092] Although the present description focuses on a digital brain using brain regions in a network, another modelling of the brain activity could be used. For instance, a 3D-map of the cortex can be used to model the brain activity of the subject, customize it to the subject and use it for training a neural network able to establish a function between a description of the frequency content of the activities of the brain and a conduction velocity.

[0093] Accordingly, the method comprises four steps, as shown in figure 1 .

[0094] In figure 1 , A illustrates that a subject-specific fractography is used to generate synthetic data based on global couplings and conduction velocities that were drawn randomly within physiological plausible ranges. C is an example of multiple global spectral features extracted from the synthetic data. D represents a machine learning used to learn an invertiblefunction that would relate velocities and global coupling to spectral features. E are posterior distributions of global coupling and conduction velocities. F and G are spectral-properties of source-reconstructed magnetoencephalography data extracted for MS patients and controls. H shows the invertible function reversed, as to estimate the most likely velocities given the observed summary statistics (i.e. spectral features).Robustness analysis:

[0095] The robustness of the model is illustrated by figure 2, which shows an estimated posterior distribution for two brain parameters from the PSD of the empirical MEG data, pooled over control and MS patient groups. In figure 2, chart A shows the global coupling strength G that shows non-significant changes (p = 0.87) and chart B shows the velocity parameter V, that significantly decreases in the control group relative to the patient group (p < 0.01 ). A reduction in the patient velocities is observed.

[0096] In figure 3, panel A, are represented observed (red) and predicted (blue) PSD of MEG data. PSDs closely resemble the empirical ones.

[0097] An example of observed and predicted MEG time-series averaged across brain regions are shown in panel B (by transforming the observation from the frequency-domain to the time-domain). The result shown in panels C and D, i.e. inferred posterior distributions for the global coupling strength G, and the velocity V, respectively, given PSD features (amplitude, median frequency, and total power), indicates that SBI accurately estimates a posterior distribution that narrows around the true parameters which generated the observed data. This confirms that SBI can accurately capture relevant PSD features of MEG data, including the peak, amplitude and total power of the PSDs.

[0098] We also checked the accuracy of the estimated posterior distributions, in that increasing the number of simulations from 100 to 20k led to progressively tighter posteriors (i.e., less uncertainty in the estimation, as shown in Fig 3, panels C and D, with the maximum a posteriori closer to the ground-truth (i.e., the parameters used to generate the data). The joint posterior distribution (i.e., on the same probability space) between parameters G and V indicates a high correlation of 0.75 (Fig. 3, panel E : joint posterior distribution between parameters G and V estimated from 20k simulations (correlation^.75). The ground truth parameters are shown in red, the high-probability parameters in yellow, the low-probability ones in blue.). Using a MAF, the training took around 3 mins, whereas generating 10000 samples from the posterior took only 1 sec. We then performed sensitivity analysis using the estimated posterior from observation, which indicated that the model is more sensitive to the parameter V than G, as conveyed by the larger eigenvalue (by around 7 orders of magnitude), which reflects a larger gradient of the posterior for the averaged velocity (panel F : sensitivityanalysis using the estimated posterior, indicating stronger model sensitivity to V than to G - the Eigenvalues for G and V are 1 .2e-05, and 8.2e-05, respectively).

[0099] In terms of computational cost, each model simulation took around 12 sec, using just-in-time compilation in Python, which can be easily run on multiple CPU / GPU cores. Using a MAF, the training took around 3 mins, whereas generating 10000 samples from the posterior took only 1 sec. We then performed sensitivity analysis using the estimated posterior from observation, which indicated that the model is more sensitive to the parameter V than G, as conveyed by the larger eigenvalue (by around 7 orders of magnitude), which reflects a larger gradient of the posterior for the averaged velocity (see figure 3, panel F).Application to a clinical cohort

[0100] The invention has applications in the medical domain, in particular for prediction or monitoring of neuro-degenerative deceases, such as multiple sclerosis.

[0101] Delays in conduction velocity of the subjects' brain are responsible for observed specific changes in power spectrum obtained from EEG or ECG.

[0102] Given that damage to the myelin in MS subjects induces delays in conduction velocity, it is hypothesized that changes in power spectrum should capture the effects of the changes in the myelin and, thence, be good predictors of clinical disability. The instant invention provides a method which makes it possible to validate this hypothesis.

[0103] In figure 4, demographic features and clinical data are combined with the resulting most likely velocities obtained from an embodiment of the method. J shows a picture of an example of the magnetic resonance of a patient with MS, showing the typical lesions. The total volume of the lesions is included, along with clinical data and demographics, in K. The outcome is a multilinear model that predicts the individual clinical disability, as measured by the EDSS.

[0104] More In detail, in a dataset made of 18 patients with MS and 20 controls, all subjects undergo both a magnetoencephalography recording and an MRI recording. Then the recordings are scanned for estimation of the tractography. This gives access to the individual connectomes and to the large-scale activity - i.e., the source-reconstructed MEG signals - for each subject.

[0105] MS was diagnosed according to the revised 2017 McDonald criteria, with exclusion criteria defined as age < 18 years, recent clinical relapse and / or steroid therapy, i.e., 3 months before the study, use of illicit drugs, stimulants, amphetamines, barbiturates, and cannabinoids; history of central nervous system disorders other than MS, severe mental illness, inability to understand and complete “patient reported outcomes” and cognitive evaluation, or inability to undergo the MRI scan.

[0106] All patients underwent a neurological clinical examination, Expanded Disability Status Scale - EDSS - scoring, the Symbol Digit Modalities Test - SDMT - to measure cognitive impairment, the Fatigue Severity Scale - FSS -, and the Beck Depression Inventory - BDI - to assess depression.

[0107] The hypothesis is validated thanks to the invention by carrying out the following tests.

[0108] An MEG recording is carried out on each participant with the same MRI scanner. The previously described MEG processing is performed.

[0109] An MRI recording is performed after the MEG recording. A processing of the diffusion MRI data is carried out as described above. A diffusion-tensor model is fitted at each voxel, and fiber tracts are generated over the whole brain by deterministic fractography using the Fiber Assignment by Continuous Tracking - FACT - algorithm described above.

[0110] To obtain the corresponding patient-specific ROI sets, each participant's FA volume is spatially normalized to the FA template provided by FSL using SPM12, and the resulting normalization matrices were inverted and applied to the ROI set.

[0111] The recordings are used to customize the digital brain model with neuronal activity data recorded on each subject.

[0112] For each subject, an invertible subject-specific function, giving the description of the frequency content of the activities of the subject's brain from an input value of global coupling parameter and average conduction velocity, is obtained by using 20000 random values and their corresponding 20000 descriptions of the frequency content of the activities.

[0113] The function is inverted into a function able to output a value of global coupling parameter and average conduction velocity.

[0114] The observed PSD on each subject is input in the continuous function. Most likely conduction velocities are thus obtained for each subject.

[0115] We then moved on to test the validity of the inferred individual delays in terms of prediction of the clinical disability. We build a multilinear model to predict the individual clinical disability (as estimated by the EDSS) based on the demographics, the clinical data and the total lesion load. Specifically, we consider the EDSS as a dependent variable, while age, education level, gender, disease duration and lesion load (i.e., the total volume of the lesions) are considered as predictors. Multicollinearity was assessed through the variance inflation factor - VIF. Figure 5 shows a multilinear regression model with leave-one-out cross validation (LOOCV) performed to test the capacity of the estimated speed to predict the EDSS scores in MS patients.

[0116] We found that the model performs well at predicting individual disability (R2 = 0.595, Adjusted R2 = 0.41 , see panels A : variance explained by the additive model including five variables (i.e., gender, age, disease duration, lesion load and estimated speed). Addingthe inferred velocities to the model significantly improves the predictive power over the EDSS (p=0.028). The predicted vs. empirical EDSS values are shown in panel B.

[0117] To test for the generalizability of our model, we used a leave-one-out cross validation (LOOCV) scheme, whereby the model is trained, at each iteration, by excluding one patient, and then the EDSS is predicted for the patient that was not included in the training data. The results of the cross-validation are shown in the lower row of Fig 5 (i.e., the average results across all iterations are reported). In this case, the average R2 is 0.61 , adj R2 = 0.40, confirming the predictive power of the model. Again, adding the estimated speed to the model significantly improved the predictive power (p=0.0417). Panels C and F show the distribution of the standardized residuals, confirming the appropriateness of using a linear model in this context. Finally, the values of Variance Inflation Factor (VIF) were always below two for all predictors, showing that multicollinearity is unlikely to affect the model.

[0118] Accordingly, the invention provides a mechanistic account of the supposed link between the modifications in the individual large-scale frequency spectrum (i.e., our observed quantity) and the conduction delays (which cannot be observed directly in MS at the wholebrain level).

[0119] The presence of such a mechanistic link allows us to use this model to infer the individual average delays, in the sense of performing Bayesian simulation-based inference in order to estimate the most likely individual delays, given the observed subject-specific large- scale frequency changes.

[0120] In figure 7, are shown the steps of a method 200 for obtaining an indicator of the progression of damage to the myelin sheath of a subject's neurons.

[0121] Step 210: comparing an earlier inferred conduction velocity of the subject's brain at an earlier time point to a later inferred conduction velocity of the subject's brain at a later time point. The earlier and later inferred conduction velocities are obtained by carrying out the method 100 of figure 6;

[0122] Step 220: providing an indicator of progression as an increasing function of the comparison between the consecutive timepoints.

[0123] In figure 8, are shown the steps of a method 300 for monitoring the efficacy of a medical treatment for multiple sclerosis.

[0124] Step 310: obtaining an indicator of the progression of damage to the myelin sheath from a recorded description of the frequency content of the activities of a subject's brain. Obtaining the indicator results from carrying out the method of figure 7;

[0125] Step 320: selection depending on the indicator of progression:

[0126] -option 322: the indicator of progression is low, meaning that the treatment is efficient or

[0127] -option 324: the indicator of progression is high, indicating a lack of efficacy of the therapy.

[0128] The invention is not limited to the disclosed embodiments.

Claims

Claims1. Computer-implemented method (100) for inferring a conduction velocity of a subject's brain from a recorded description of the frequency content of the activities of the subject's brain, characterized in that it comprises the steps of:(i) customizing (1 10) a digital brain model with neuronal activity data recorded on the subject, the digital brain model being capable of simulating realistic personalized activities of the subject's brain in the form of time-series,(ii) inputting (120) n plausible values of average conduction velocity into the customized digital brain model to obtain time-series in order to characterize n descriptions of the frequency content of the activities of the simulated brain,(iii) establishing (130) a backward function from the calculated descriptions to the inputted values,(iv) applying (140) the backward function to the recorded description of the frequency content of the activities of the subject's brain to obtain the average conduction velocity of the subject's brain.

2. Method according to claim 1 , wherein the digital brain model is a network of a plurality of brain regions, each brain region being modelized by a local area model, the local area models being connected to each other in the network configured according to a fractography carried out on the subject.

3. Method according to claim 2, wherein any inputted average conduction velocity value is combined with any global coupling value.

4. Method according to claim 1 , wherein the digital brain model is based on a 3D-map of the cerebral cortex.

5. Method according to any one of claims 1 , 2, 3 and 4, wherein the backward function is established using an artificial neural network.

6. Method according to any one of claims 1 , 2, 3, 4 and 5, wherein the backward function is obtained by first establishing an invertible forward function from the inputted values to the calculated descriptions, then inverting the forward function from the descriptions to the conduction velocities.

7. Method according to any one of claims 1 , 2, 3, 4, 5 and 6, wherein the description of the frequency content of the activities of a brain is a power spectra description.

8. Method according to any one of claims 1 , 2, 3, 4, 5, 6 and 7, wherein the neuronal activity data contains non-invasive magneto-encephalography data or electro-encephalography data.

9. Method according to any one of claims 1 , 2, 3, 4, 5, 6, 7 and 8, wherein in step (120), 15,000 to 25,000 plausible values, preferably 20,000, are inputted into the customized digital model.

10. Method according to any one of claims 1 , 2, 3, 4, 5, 6, 7, 8 and 9, wherein the plausible values are randomly set in a range of 0.5ms to 20ms.

11. Computer program product, characterized in that it comprises instructions which, when the program is executed by a computer, cause the computer to perform the method according to any one of claims 1 to 10.

12. Data processing system, characterized in that it comprises a processor configured to perform the method (100) according to any one of claims 1 to 10.

13. Computer-readable storage medium, characterized in that it comprises instructions which, when executed by a computer, cause the computer to perform the method (100) according to any one of claims 1 to 10.

14. Method for obtaining an indicator (200) of the progression of damage to the myelin sheath of a subject's neurons, comprising the steps of:- comparing (210) an earlier inferred conduction velocity of the subject's brain at an earlier time point to a later inferred conduction velocity of the subject's brain at a later time point;- providing (220) an indicator of progression as an increasing function of the comparison between the consecutive timepoints, characterized in that the earlier and later inferred conduction velocities are obtained by carrying out the method (100) according to any one of claims 1 to 10.

15. Method for monitoring (300) the efficacy of a medical treatment for multiple sclerosis, comprising the steps of:- obtaining (310) an indicator of the progression of damage to the myelin sheath from a recorded description of the frequency content of the activities of a subject's brain,- indicating that the treatment is efficient when the indicator of progression is low (322) or indicating a lack of efficacy of the therapy when the indicator of progression is high (324), characterized in that obtaining the indicator results from carrying out the method (200) according to claim 14.

Images