Méthodes quantitatives en neurosciences

Cours NSC-2006, année 2015
Laboratoire d'analyse de données multidimensionnelle-Réponses
*Pierre Bellec, Yassine Ben Haj Ali*

Objectif:

Ce laboratoire a pour but de vous initier à la manipulation d’informations multidimensionnelles avec Matlab. Nous allons pour cela analyser des données électrophysiologiques de décharges d’activité neuronale. Nous allons effectuer différentes opérations visant à visualiser, résumer et modéliser ces données. Les données sont tirée de l’expérience de Georgopoulos1982 sur l’encodage neuronale du mouvement du bras chez un macaque avec des implants neuronaux. L’animal commence l’expérience en fixant un curseur au centre d’une cible, ensuite il doit rejoindre des cibles périphériques apparaissent dans une des 8 directions arrangé en cercle. Une fois la cible apparue, l’animal doit attendre ( 100-1500 ms) le signal de départ avant d’aller rejoindre la cible pour une durée de 500ms, ensuite il retourne au point de départ (au centre). Cette séquence de mouvement est appelée essai et dans cette expérience il y en a 47. Le but de l’expérience de Georgopoulos et collègues est de déterminer l’orientation spatiale préférentielle du neurone en question dans la région MI, et qu’il est possible de prédire la direction du mouvement à partir d’enregistrements physiologiques. Leurs résultats indiquent qu’il y a bel et bien une préférence vers les angles de mouvement entre 90 et 180 degrés. Durant ce travail nous allons reproduire certaines des analyses de données et la visualisation des résultats de cette expérience.

Pour réaliser ce laboratoire, il est nécessaire de récupérer les resources suivantes sur studium:

Chap17_Data.mat: le jeu de données tiré de Georgopoulos1982.
Les scripts diagramme_dispersion.m et diagramme_dispersion_essais.m pour la Section 1.
Les scripts histogramme_essai1.m, histogramme_essais.m pour la Section 2.

Notez que le laboratoire est noté. Il faudra remettre un rapport détaillé incluant une réponse pour l’ensemble des questions numérotées ci dessous. Chaque réponse fera typiquement quelques lignes, incluant du code et une figure si demandé dans l’énoncé.

Ne pas tenir compte et ne pas executer cette partie du code:



In [39]:

    
%matplotlib inline
from pymatbridge import Matlab
mlab = Matlab()
mlab.start()
%load_ext pymatbridge









    



Starting MATLAB on ZMQ socket ipc:///tmp/pymatbridge
Send 'exit' command to kill the server
.MATLAB started and connected!
The pymatbridge extension is already loaded. To reload it, use:
  %reload_ext pymatbridge

Section 1 : Diagramme de dispersion

Nous allons commencer par effectuer un diagramme de dispersion (scatter plot) de l’activation d’un neurone tout au long de la durée d’un essai. Voir le script pour suivre ces étapes.

Commençons par charger les données:



In [40]:

    
%%matlab
load('Chap17_Data')

La commande `whos` nous permet de déterminer quelles variables sont disponibles dans l’espace de travail, ainsi que leur type.



In [41]:

    
%%matlab
whos









    





  Name               Size                Bytes  Class     Attributes

  ans                1x1                     8  double              
  direction        158x1                  1264  double              
  ftheta             1x1                  1266  inline              
  go               158x1                  1264  double              
  instruction      158x1                  1264  double              
  spike              1x47                29624  struct              
  theta_chap         1x3                    24  double              
  unit               1x143            11872990  struct              
  x                  1x301                2408  double              
  y                  1x301                2408  double

Quelles variables sont présentes? Quel est le type de la variable spike? Quelle est sa taille?

Les variables sont `direction`, `go`, `intruction`, `spike` et `unit`. La variable `spike` est de type structure et de taille 1x47

La variable spike contient les temps des potentiels d’action détectés pour un neurone. Chaque entrée de la structure contient les données d’un essai différent. Il est possible de lister les champs de la structure avec la commande fieldnames:



In [42]:

    
%%matlab
fieldnames(spike)









    





ans = 

    'times'

Le champ "spikes(1).times" contient les temps de décharges de potentiels d’action pour
le premier essai. La commande permet de déterminer la taille de ce
vecteur, c’est à dire le nombre de décharges détectées:



In [43]:

    
%%matlab
size(spike(1).times)

Combien y-a-t-il eu de décharges pour l’essai 2? pour l’essai 10?



In [44]:

    
%%matlab
size(spike(2).times) %nb de décharges pour l'essai 2

Il y a 51 décharges pour l'essai 2.



In [45]:

    
%%matlab
size(spike(10).times) %nb de décharges pour l'essai 10

Il y a 85 décharges pour l'essai 10.

La commande suivante va présenter l’ensemble des temps de décharge pour l’essai 1.

        >> spike(1).times
           -0.9893
           -0.9402
           -0.9158
        (...)

Quelle est l’unité probable de ces temps? Pourquoi y-a-t-il des valeurs négatives?

L'unité est en secondes et les valeurs négatives indiquent un temps de décharge qui précède le signal de départ le de l'expérience

On extrait les temps de décharges des deux premiers essais dans deux variables `t1` et `t2`:



In [46]:

    
%%matlab
 t1 = spike(1).times; 
 t2 = spike(2).times;

On ouvre une nouvelle fenêtre, dédiée à la visualisation:



In [47]:

    
%%matlab
figure 
hold on

Attention à la deuxième instruction! Elle permet de dessiner
plusieurs objets sur une même figure, l’un à la suite de l’autre,
sans ré-initialiser la figure.

Maintenant on va tracer la première ligne du diagramme. Notez que le nombre de décharges dans l’essai 1 est . On applique une boucle :



In [48]:

    
%%matlab
for num_temps = 1:length(t1) 
    line([t1(num_temps) t1(num_temps)], [0 1]) 
end

Notez l’utilisation de la commande line.

On va maintenant ajouter un label sur l’axe des x et y .

Sauvegardez la figure dans un fichier png (utilisez la commande print), sous le nom figure_dispersion.png



In [49]:

    
%%matlab
for num_temps = 1:length(t1) 
    line([t1(num_temps) t1(num_temps)], [0 1]) 
end
xlabel('Temp (sec)'); 
%Idem pour l’axe des y:
ylabel('Essai #')
%Enfin, on fixe les limites de l’axe des y
ylim([0 3])
% save the result
print('figure_dispersion.png','-dpng')

Faire une nouvelle figure où chaque barre du diagramme a une hauteur de 0.5, plutôt que 1. Sauvegardez ce fichier sous le nom "figure_dispersion_lignes.png".



In [50]:

    
%%matlab
for num_temps = 1:length(t1) 
    line([t1(num_temps) t1(num_temps)], [0 0.5]) 
end
xlabel('Temp (sec)')
%Idem pour l’axe des y:
ylabel('Essai #')
%Enfin, on fixe les limites de l’axe des y
ylim([0 3])

Vous allez compléter, à partir du fichier , les 4 lignes de code manquante à l’interieure de la boucle pour tracer tous les essais (47) dans une même figure. Le résultat ressemblerait à la figure suivante:



In [51]:

    
%%matlab
% Charger les donnees
load('Chap17_Data')
% Preparer une figure
figure 
% permettre la superposition de plusieurs graphiques dans la meme figure
hold on 
% Donner un label à l'axe des x
xlabel('Temp (sec)'); 
% Donner un label à l'axe des y
ylabel('Essai #');
% Ajuster les limites de l'axe des y
ylim([0 length(spike)]);

for num_spike = 1:length(spike) %faire une boucle pour tout les essaies
    t = spike(num_spike).times; %definir la variable pour chaque essai
    for num_temps=1:length(t) %faire une boucle pour tous les points temps
        line([t(num_temps) t(num_temps)], [0+(num_spike-1) 1+(num_spike-1)]); %dessiner une line pour chaque point temps t1(i) avec longueur de 1
    end
end

Section 2 : Histogramme

Nous allons continuer l’exploration des données à travers un histogramme qui résumerait la somme des activations dans un intervalle de temps donné. Voir le script histograme_essai1.m pour reproduire les commandes suivantes:

Commençons par nettoyer l’espace de travail:



In [52]:

    
%%matlab
clear

Chargeons de nouveau les données:



In [53]:

    
%%matlab
load('Chap17_Data')

Définissons les bords et le pas des catégories de l’histogramme



In [54]:

    
%%matlab
centres = [-0.95:0.1:0.95];

Initialiser une matrice de zéros dont la longueur est égale au nombre d’intervalles:



In [55]:

    
%%matlab
histo = zeros(1,length(centres));

Récupérez le nombre de décharges par intervalle de temps dans l’essai numéro 1, à l’aide avec la fonction .



In [56]:

    
%%matlab
histo = hist(spike(1).times,centres);

Examinez le contenu de la variable `histo`. Quelle est sa taille ? Son minimum, maximum, sa moyenne (voir les fonctions matlab `min, max, mean` ).



In [57]:

    
%%matlab
whos histo % elle est de taille 1x20









    





  Name       Size            Bytes  Class     Attributes

  histo      1x20              160  double



In [58]:

    
%%matlab
min(histo)
max(histo)
mean(histo)









    





ans =

     1


ans =

     6


ans =

    2.6000

Dessinez l’histogramme avec la fonction `bar`.



In [59]:

    
%%matlab
bar(centres,histo);
 %Ajuster les limites de l’axe des x
xlim([-1.1 1]);
xlabel('Temps (sec)');  %Donner un label à l’axe des x
ylabel('# essai');%Donner un label à l’axe des y

Reprenez le code du fichier et remplissez la boucle afin de réaliser un histogramme pour l’ensemble des essais



In [60]:

    
%%matlab
%Charger les donnees
load('Chap17_Data')

% Definir les centres des intervalles pour l'histogramme
centres = [-0.95:0.1:0.95];

% Initialiser une matrice de zéros histo dont la longueur est égale au nombre d'intervalles:
histo = zeros(length(centres),1);

% Faire une boucle à travers tous les essais et recuperer le nombre de decharges par intervalle avec la fonction histc 
for jj = 1:47
    histo=histo+histc(spike(jj).times,centres);
end

% Dessiner l'histograme avec la fonction bar
bar(centres,histo);

%Ajuster les limites de l'axe des x
xlim([-1.1 1]);

%Donner un label à l’axe des x
xlabel('Temps (sec)');  

%Donner un label à l’axe des y
ylabel('# essai');

Section 3 : Régression

Nous allons maintenant implémenter une régression à l’aide de Matlab.

Commençons par nettoyer l’espace de travail:



In [61]:

    
%%matlab
clear

Maintenant nous allons récuperer les données de “notes” du cours.



In [62]:

    
%%matlab
x = [ 165 165 157 170 175 165 182 178 ]';
y = [  47  56  49  60  82  52  78  90 ]';

Quels est la taille et le contenu des vecteurs x et y? A quoi sert l'opération `'` ?



In [63]:

    
%%matlab
whos









    





  Name      Size            Bytes  Class     Attributes

  x         8x1                64  double              
  y         8x1                64  double

Les variables x et y sont tout deux des vecteurs de 1 colonne et 8 lignes. L'operation `'` transposer un vecteur ligne en colonne ou le contraire.

On définit une nouvelle fonction en ligne:



In [64]:

    
%%matlab
ftheta = inline('theta(1)+theta(2)*x','theta','x');

Quel est le type de la variable `ftheta`?



In [65]:

    
%%matlab
whos ftheta









    





  Name        Size            Bytes  Class     Attributes

  ftheta      1x1              1238  inline

`ftheta` est une fonction en ligne (*inline*).

Estimez les coefficients de régression à l’aide de la fonction :



In [66]:

    
%%matlab
theta_chap = nlinfit(x, y, ftheta, [1 1] );

Quelles sont les valeurs de `theta_chap`? A quoi sert l’argument `[1 1]`? Essayez de reproduire l’estimation avec d’autres valeurs pour cet argument, est-ce que cela affecte le résultat?

Le paramètre `theta_chap` vaut:



In [67]:

    
%%matlab
theta_chap









    





theta_chap =

 -237.5729    1.7794

L'argument `[1 1]` est une valeur initiale de la méthode qui cherche la valeur `theta_chap`. En répétant l'expérience pour plusieurs valeurs (`[2 2], [30 30], [-30 -30]`) on voit que la résultat `theta_chap` ne semble pas dépendre ici de ce paramètre.

Maintenant représenter le résultat de la régression.



In [68]:

    
%%matlab
figure
plot(x,y,'b.');
hold on
plot(x,ftheta(theta_chap,x),'r');

Utilisez la fonction `ylim` pour changer les limites de l’axe y de 40 à 95. Ajouter le label `taille` sur l’axe des x avec la commande `xlabel`, et le label `poids` sur l’axe des y avec la commande `ylabel`. Faites une sauvegarde de cette image, dans un fichier `regression_notes.png`.



In [69]:

    
%%matlab
figure
plot(x,y,'b.');
hold on
plot(x,ftheta(theta_chap,x),'r');
ylim([40 95])
xlabel('taille')
ylabel('poids')
print('regression_notes.png','-dpng')

Maintenant nous allons ajuster une courbe plus complexe, un cosinus. On commence par simuler des données:



In [70]:

    
%%matlab
clear
x = 0:0.1:30;
y = cos(x) + randn(1,301);

Quelle est la taille de x? La taille de y?



In [71]:

    
%%matlab
size(x)
size(y)

Les variables `x` et `y` sont des vecteurs lignes de longueur 301.

A quoi sert la fonction `randn` (utilisez la commande `help` ).



In [72]:

    
%%matlab
help randn









    





 RANDN Normally distributed pseudorandom numbers.
    R = RANDN(N) returns an N-by-N matrix containing pseudorandom values drawn
    from the standard normal distribution.  RANDN(M,N) or RANDN([M,N]) returns
    an M-by-N matrix. RANDN(M,N,P,...) or RANDN([M,N,P,...]) returns an
    M-by-N-by-P-by-... array. RANDN returns a scalar.  RANDN(SIZE(A)) returns
    an array the same size as A.
 
    Note: The size inputs M, N, P, ... should be nonnegative integers.
    Negative integers are treated as 0.
 
    R = RANDN(..., 'double') or R = RANDN(..., 'single') returns an array of
    normal values of the specified class.
 
    The sequence of numbers produced by RANDN is determined by the settings of
    the uniform random number generator that underlies RAND, RANDN, and RANDI.
    RANDN uses one or more uniform random values to create each normal random
    value.  Control that shared random number generator using RNG.
  
    Examples:
 
       Example 1: Generate values from a normal distribution with mean 1
        and standard deviation 2.
          r = 1 + 2.*randn(100,1);
 
       Example 2: Generate values from a bivariate normal distribution with
       specified mean vector and covariance matrix.
          mu = [1 2];
          Sigma = [1 .5; .5 2]; R = chol(Sigma);
          z = repmat(mu,100,1) + randn(100,2)*R;
 
       Example 3: Reset the random number generator used by RAND, RANDI, and
       RANDN to its default startup settings, so that RANDN produces the same
       random numbers as if you restarted MATLAB.
          rng('default');
          randn(1,5)
 
       Example 4: Save the settings for the random number generator used by
       RAND, RANDI, and RANDN, generate 5 values from RANDN, restore the
       settings, and repeat those values.
          s = rng
          z1 = randn(1,5)
          rng(s);
          z2 = randn(1,5) % z2 contains exactly the same values as z1
 
       Example 5: Reinitialize the random number generator used by RAND,
       RANDI, and RANDN with a seed based on the current time.  RANDN will
       return different values each time you do this.  NOTE: It is usually
       not necessary to do this more than once per MATLAB session.
          rng('shuffle');
          randn(1,5)
 
    See <a href="matlab:helpview([docroot '\techdoc\math\math.map'],'update_random_number_generator')">Updating Your Random Number Generator Syntax</a> to use RNG to replace
    RANDN with the 'seed' or 'state' inputs.
 
    See also RAND, RANDI, RNG, RANDSTREAM, RANDSTREAM/RANDN

    Overloaded methods:
       RandStream/randn

    Reference page in Help browser
       doc randn

La commande `randn` permet de simuler du bruit suivant une distribution normale (Gaussienne) de moyenne nulle et de variance 1.

Générer un graphe de la relation entre x et y, et sauvegardez cette image dans un fichier `donnees_cosinus.png`.



In [73]:

    
%%matlab
figure
plot(x,y,'.')
print('donnees_cosinus.png','-dpng')

On va maintenant définir une fonction de trois paramètres:



In [74]:

    
%%matlab
ftheta = inline('theta(1)+theta(2)*cos(x-theta(3))','theta','x');

Quelle est la valeur de la fonction `ftheta` pour `theta=[0 1 1]` et `x=0` ?



In [75]:

    
%%matlab
ftheta([0 1 1],0)

Estimez les coefficients de régression à l’aide de la fonction :

Quelles sont les valeurs de `theta_chap` ?



In [76]:

    
%%matlab
theta_chap = nlinfit(x, y, ftheta, [0 1 1] )









    





theta_chap =

    0.0088    1.1351   -0.0188

A quoi sert l’argument `[0 1 1]` ? Essayez de reproduire l’estimation avec d’autres valeurs pour cet argument, est-ce que cela affecte le résultat?

L'argument [0 1 1] jour le même rôle que le paramètre [1 1] à la question 4. En utilisant une valeur différente du paramètre (par exemple [0 1 10]) on trouve une valeur différente pour `theta_chap`. C'est parce que la fonction `ftheta` est périodique, et il existe donc une infinité de valeurs d'entrée qui donnent la même sortie.

Maintenant représenter le résultat de la régression.



In [77]:

    
%%matlab
figure
plot(x,y,'b');
hold on
plot(x,ftheta(theta_chap,x),'r');

Faites une sauvegarde de cette image, dans un fichier .