Séance 1 : Modéliser des flux protéiques dans la cellule grâce à l'algèbre linéaire

import numpy as np # numpy est la bibliothèque de calcul numérique en python. 

import matplotlib.pyplot as plt # Matplotlib permet de tracer des graphes.

# la commande jupyter ci-dessous permet d'afficher les graphes de matplotlib directement dans le notebook.
%matplotlib inline

from matplotlib.patches import Polygon # objet permettant de tracer des polygones

plt.rc('text', usetex=True) # Configure matplotlib pour utiliser TeX pour le rendu des textes 

Une histoire de flux

Le transport des protéines dans la cellule eucaryote est un phénomène complexe faisant appel à de nombreux compartiments (reticulum endoplasmique, noyau, appareil de Golgi, vésicules, mitochondries, plastes, membrane) et à des flux de matière très controlés qui les relient.

On s'intéresse au suivi de protéines fluorescentes marquées à la GFP dans la cellule. On souhaite décrire l'évolution de leur localisation dans les différents compartiments cellulaires, sachant que ces protéines marquées peuvent être transportées entre compartiments, ou dégradées, et que de nouvelles protéines peuvent être synthétisées.

Pour cela, nous proposons de réaliser un modèle linéaire en compartiments.

Exercice 0:

Commençons en une dimension par considérer l'évolution de la quantité totale de protéines dans la cellule, $u(t)$. On suppose qu'à chaque pas de temps, cette quantité est multipliée par une constante $a$ (système linéaire - suite géométrique), tel que $$u(t+1)= a u(t)$$

Représenter la trajectoire pour $u(0)=1$ et a=1, a=1.15 ou a=0.95. Interpréter. (Indice : utiliser plt.plot).

# Correction exercice 0 :

def display_traj_suite(u0=0.5,T=20,a=1):
    u=[u0,]
    for i in range(T):
        u.append(a*u[-1])

    plt.plot(u, 'o')

    
display_traj_suite(u0=1,T=15,a=1.15)
display_traj_suite(u0=1,T=15,a=1)
display_traj_suite(u0=1,T=15,a=0.95)

plt.legend(["1.15","1","0.95" ])

<matplotlib.legend.Legend at 0x119fd1d30>

Rappels sur les matrices

Une matrice $A$ de taille $m \times n$ à coefficient dans $\mathbb K$ est un ensemble d'éléments de $\mathbb K$ ordonnés en $m$ lignes et $n$ colonnes.On note $a_{ij}$ l'élément à la $i\text{-ème}$ ligne et $k\text{-ème}$ colonne.

\begin{equation} A = (a_{ij})_{1\leq i\leq m; 1 \leq j \leq n} \end{equation}

On peut aussi commencer l'indexation à 0 pour coller à la représentation informatique: $A = (a_{ij})_{0 \leq i < m; 0 \leq j < n}$

Exemple en taille $3 \times 3$:

\begin{equation} A = \left [ \begin{array}{ccc} a_{00} & a_{01} & a_{02} \\ a_{10} & a_{11} & a_{12} \\ a_{20} & a_{21} & a_{22} \\ \end{array} \right ] \end{equation}

# Les matrices peuvent être représentées par des objets de type numpy.array. 
# Créer un array peut se faire à  partir d'une liste de liste:
A = np.array([[10,6,6],
              [8,2,4],
              [8,4,2]])

# Accéder à  $a_{ij}$ se fait par indexation:
A[0,1] # => 6

# Il existe des constructeurs pour des matrices particulières:
I = np.identity(3) # Matrice identité de taille 3
N = np.zeros((3,3)) # Matrice de zeros de taille 3
X,Y = np.meshgrid([1,2,3],[1,2,3]) # Matrices donnant le produit cartésien de deux vecteurs (permet d'obtenir un maillage).

x = np.array([1,2,3])
# Le produit matriciel pour les np.array peut se faire avec la méthode .dot:
A.dot(x) # Seule possibilité pour Python < 3.5
# Ou l'opérateur @:
A @ x #Python >= 3.5
assert np.all(A.dot(x) ==  A @ x)

# Le produit matriciel n'est défini que pour des matrice de dimensions compatibles. 
# Exemple, le produit: A * B necessite que A ait autant de colonnes que B ait de lignes. 
# Le résultat est une matrice carrée. 
# (n,m) * (m,k) -> (n,k)

# Par défaut les array numpy de dimension 1 sont orientés à  chaque opération pour que le produit matriciel fonctionne :
x = np.array([1,2,3])
print("value: {}, dim: {}, shape: {}, A@x: {}, x@A {}".format(x, x.ndim, x.shape, A@x, x@A))

# Ainsi, Pour deux arrays de dimension 1, @ est le produit scalaire.
x@x # équivalent à  <x,x> ou x^T * x si x est un vecteur colonne.

value: [1 2 3], dim: 1, shape: (3,), A@x: [40 24 22], x@A [50 22 20]

14

# La syntaxe np.newaxis permet d'augmenter sa dimension. Maintenant, l'orientation ligne ou colonne est importante.
b = x[:,np.newaxis] # `b` est un vecteur colonne (3 lignes, 1 colonne)
print("value: {}, dim: {}, shape: {}, A@b: {}".format(b, b.ndim, b.shape, A@b))

c = x[np.newaxis,:] # `c` est un vecteur ligne (1 ligne, 3 colonnes)
print("value: {}, dim: {}, shape: {}, b@A:{}".format(c, c.ndim, c.shape, c@A))

c@b #=> un scalaire (produit scalaire des deux vecteurs <b,c>).
b@c #=> une matrice 3x3. 

value: [[1]
 [2]
 [3]], dim: 2, shape: (3, 1), A@b: [[40]
 [24]
 [22]]
value: [[1 2 3]], dim: 2, shape: (1, 3), b@A:[[50 22 20]]

array([[1, 2, 3],
       [2, 4, 6],
       [3, 6, 9]])

# NB : l'opérateur "*" avant un vecteur permet de "déballer" ses coordonées (nécessaire lorsque l'on veut plotter des vecteurs). 
x = [1.0, 0.1, 0.5]

print(*x)

1.0 0.1 0.5

Partie 1 : Simuler des trajectoires

On considère le modèle simplifié suivant, où les transports sont représentés en vert, les synthèses en jaune et les dégradations en rouge :

Exercice 1 :

Écrire la matrice A qui décrit ce système dynamique.

\begin{equation} \mathbf s = \left [ \begin{array}{c} s\\0\\0 \end{array} \right ] , \; \mathbf d = \left [ \begin{array}{c} 0\\0\\d\end{array} \right ] , \; \phi = \left [ \begin{array}{ccc} 0 & r & 0 \\ f & 0 & g \\ m & e & 0 \\ \end{array} \right ] \end{equation}\begin{equation} \mathbf A = \mathbf I + \mathbf I \mathbf s - \mathbf I \mathbf d + \phi \end{equation}\begin{equation} \mathbf A = \left [ \begin{array}{ccc} 1-f-m+s & r & 0 \\ f & 1-e-r & g \\ m & e & 1-g-d \\ \end{array} \right ] \end{equation}

Exercice 2 :

1) Écrire une fonction qui construit la matrice A tout en vérifiant la conservation de la matière. (Indice : utiliser assert).

2) Simuler la trajectoire de ce modèle simplifié, pour f=0.4, m=0.2, r=0.5, e=0.4, g=0.1, s=0.5, et d=0. Tracer le graphe du contenu de chaque compartiment au cours du temps, sachant qu'au départ, la protéine est seulement présente en concentration 1 dans le Golgi. (Indice : utiliser plt.subplots).

3) Représenter différents types de trajectoires, en faisant varier le ratio synthèse/dégradation. (Indice : utiliser plt.subplots).

# Correction exercice 2 - question 1 :

def transport_3comp(f=.4,m=.2,r=.5,e=.4,g=.1,s=0, d=0):
    """ Matrice de la dynamique du modèle à  3 compartiments.
    Args: 
        f(float): flux golgi->endo.
        m(float): flux golgi->surface.
        r(float): flux endo->golgi.
        e(float): flux endo->surface.
        s(float): synthese dans le golgi.
        d(float): degradation a la surface
    Returns: un np.array de taille 3x3. 
    """
    # Programmation defensive ! 
    # Les flux sortants doivent être <100% 
    # pour assurer la conservation de la matière. 
    assert f+m <= 1, 'flux de sortie du golgi >100%'
    assert e+r <= 1, 'flux de sortie des endosomes >100%'
    assert g+d <= 1, 'flux de sortie de la surface cellulaire >100%'
    A = np.array([[1-f-m+s,r,0],
                    [f,1-e-r,g],
                    [m,e,1-g-d]])
    return A


# Matrice de la dynamique
A = transport_3comp(s=0.5)
A

array([[0.9, 0.5, 0. ],
       [0.4, 0.1, 0.1],
       [0.2, 0.4, 0.9]])

# Correction exercice 2 - question 2 :

def simulate(A, u0, T=10):
    """Simule le modèle u_n = A^n u_0
    Args:
        A (np.array): Matrice de la dynamique
        u0 (iterable): Condition initiale
        T (int): Nombre de pas de temps
    Renvoie la trajectoire sous forme de np.array: traj[n,i] = (u_n)_i
    """
    traj = [u0, ]
    for t in range(T):
        traj.append(A@traj[-1])
    return np.array(traj)

def display_traj3comp(traj, ax=None, labels=['Golgi','Endosomes', 'Surface']):
    """Affiche la trajectoire du modèle de transport à 3 compartiment"""
    if ax is None:
        fig, ax = plt.subplots(figsize=(12,4))
    ax.plot(traj) # trace les lignes
    ax.legend(labels) # ajoute légendes
    for t in traj.T:
        ax.scatter(np.arange(len(t)),t)  # superpose points


# Conditions initiales
u0 = [1,0,0]

# Temps final
T = 30

# Simulation.
traj = simulate(A, u0, T)

# Affichage
display_traj3comp(traj)

# Correction exercice 2 - question 3 :

# Créer une liste avec tous les jeux de paramètres à  tester.
parametres = [("Production",{'s':0.5}), 
              ("Degradation",{'d':0.5}),
              ("Production $>$ Dégradation",{"s":2, 'd':0.5})]

# Créer la figure
fig,axes = plt.subplots(1,3, figsize=(15,4))

# Itérer sur les jeux de paramètres et les subplots (objets ax).
for (title,p),ax in zip(parametres, axes):
    A = transport_3comp(**p)
    traj = simulate(A, u0, T)
    display_traj3comp(traj, ax=ax)
    ax.set(title=title)

    
parametres = [("Dégradation $>$ Production",{"s":0.1, 'd':.8}),
              ("Dégradation $>$ Production",{"s":0.1, 'd':.8,'m':0,'g':0,'r':0,'f':1,'e':1}),
              ("Périodique",{'m':1,'g':1,'r':1,'f':0,'e':0}),]

# Créer la figure
fig,axes = plt.subplots(1,3, figsize=(15,4))

# Itérer sur les jeux de paramètres et les subplots (objets ax).
for (title,p),ax in zip(parametres, axes):
    A = transport_3comp(**p)
    traj = simulate(A, u0, T)
    display_traj3comp(traj, ax=ax)
    ax.set(title=title)
    

    
plt.savefig('trajectoires.pdf')

Partie 2 : matrices $2\times 2$, des endomorphismes sur le plan

Pour encore simplifier et essayer de se construire une intuition géométrique on va passer à des matrices $2 \times 2$.

\begin{equation} A = \left [ \begin{array}{cc} a_{00} & a_{01} \\ a_{10} & a_{11} \\ \end{array} \right ] \end{equation}

Représenter l'effet d'une matrice

Comment représenter une application en fonction de ses espaces de départ et d'arrivée

Espace de départ	Espace d'arrivée	Représentation graphique	Matplotlib
$\mathbb R$	$\mathbb R$	courbe, points	plt.plot, plt.scatter
$\mathbb R^2$	$\mathbb R$	contour, surface	plt.contour, plt.contourf, plt.imshow, axes3D.plot_surface
$\mathbb R^2$	$\mathbb R^2$	champ de vecteur, écoulement	plt.quiver, plt.streamplot

Exercice 3 :

Soit la matrice \begin{equation} A = \left [ \begin{array}{cc} 1.2 & 0.46 \\ 0.45 & 0.8 \\ \end{array} \right ] \end{equation}

1) Représenter l'effet du produit matriciel sur différents vecteurs du plan. (Indice: utiliser les fonctions plt.scatter et plt.arrow. Options : plt.set, plt.legend et plt.grid).

2) Illustrer la transformation en utilisant un champs de vecteur. (Indice: utiliser les fonctions np.linspace, np.meshgrid, np.vectorize et plt.quiver).

3) Représenter l'effet du produit matriciel sur le carré unité (i.e. carré formé à partir des vecteurs de la base canonique). (Indice: utiliser la fonction plt.add_patch(Polygon())). Superposer un cercle unité (sachant que son équation correspond à $x^2+y^2=1$).

4) Interpréter ces transformations en terme d'évolution du système Réticulum-Golgi.

# Correction exercice 3 - question 1 :

A = np.array([[1.2,0.46],
              [0.45,0.8]])
x = [0.7,0.3]

# Voyons, graphiquement, l'effet d'un produit matriciel sur un point:
fig, ax = plt.subplots(figsize=(8,8))
ax.scatter(*x, label=r'$\vec{x}$')    # Ajoute un point 
ax.scatter(*A@x, label=r'$A \circ \vec{x}$')
ax.scatter(0,0, label=r'0',color='k')
ax.arrow(0,0,*x,color='C0')     # Trace une ligne entre 2 points
ax.arrow(0,0,*A@x,color='C1')

ax.set(xlim=(0,1.3*np.max(A@x)), ylim=(0,1.3*np.max(A@x))) # Limites des axes
ax.legend() # Affiche la légende
ax.grid(1)  # Quadrillage 
 
plt.savefig('produit_matriciel.pdf')

# Correction exercice 3 - question 2 :

lspace = np.linspace(1,10,5)
fig, ax = plt.subplots(figsize=(10,10))

# Meshgrid donne le produit cartésien des deux vecteurs.
# Et donc de la grille sur laquelle on veut appliquer A.
X,Y = np.meshgrid(lspace,lspace)

# np.vectorize permet de donner à une fonction quelconque (ici le lambda)
# des capacités de broadcasting (application à tous les éléments).
# Par défaut vectorize attend une fonction qui prend des scalaires et retourne des scalaires.
# Si ce n'est pas le cas il faut préciser la signature. Ici la fonction prend deux scalaires "(),()" et
# retourne un array de taille 2 "->(2)"
f = np.vectorize(lambda x,y: A@[x,y], signature='(),()->(2)')
U = f(X,Y) # On applique la fonction à tous les éléments de X,Y ! 

ax.scatter(X, Y, color='C0', label=r'$\vec{x}$')
ax.scatter(U[:,:,0],  U[:,:,1], color='C1', label=r'$A \circ \vec{x}$')
ax.quiver(X, Y, U[:,:,0]-X, U[:,:,1]-Y, angles='xy', units='xy', scale=1, color='C2', alpha=0.4)
ax.set(xlim=(0,np.max(U)),ylim=(0,np.max(U)))
ax.legend()
ax.grid()

plt.savefig('champs_vecteurs.pdf')

L'aire de l'image du carré unité est égal à la valeur absolue du déterminant. Le déterminant est négatif si l'orientation du carré unité a changé (les vecteurs de la base canonique se sont croisés).

# Correction exercice 3 - question 3 :

def display_transform(A, ax=None, title=None, name='A'):
    """Affiche l'effet de la matrice A sur le carré unité"""
    if ax is None:
        fig, ax = plt.subplots(figsize=(6,6))
    if title is not None:
        ax.set(title=title)
        
    # unit square
    square = [(0,0),(0,1),(1,1),(1,0)]

    # Plot the canoncial basis and its image 
    arrow ={'linewidth':2, 'head_width':0.05, 'length_includes_head':True}
    ax.arrow(0, 0, 0, 1, color='k', label='e', **arrow)
    ax.arrow(0, 0, 1, 0, color='k', **arrow)
    ax.arrow(0, 0, *A@np.array([0,1]), color='C0' , **arrow)
    ax.arrow(0, 0, *A@np.array([1,0]), color='C0', label=name+'e', **arrow)
        
    # Plot the unit square and its image
    Asquare = [A@x for x in square]
    p1=ax.add_patch(Polygon(square, facecolor='k', alpha=0.3))
    p2=ax.add_patch(Polygon(Asquare, facecolor='C0', alpha=0.3))
    
    ## Plot a unit circle
    u = np.linspace(-1,1,1000)
    ax.plot(u, np.sqrt(1-u**2),color='k',alpha=0.1)
    ax.plot(u, -np.sqrt(1-u**2),color='k',alpha=0.1)
    
    # Plot a grid 
    ax.grid(1)
    
    # Informative legend
    tex =  ('${}'.format(name)
            +r' = \left [ \begin{array}{cc}'
            + r"{0[0][0]:2} & {0[0][1]:2} \\ {0[1][0]:2} & {0[1][1]:2}".format(A) 
            + r'  \end{array} \right ] $')
    ax.legend((p1,p2),('$e$','${} \circ e$'.format(name)), title=tex)
    
    # Set axis limits
    ax.set(xlim=(np.floor(np.min([np.min(Asquare),-1])),
                 np.ceil(np.max([np.max(Asquare),1]))), 
           ylim=(np.floor(np.min([np.min(Asquare),-1])),
                 np.ceil(np.max([np.max(Asquare),1]))))
    return ax

display_transform(A);

plt.savefig('carré_unité.pdf')

Produit Matriciel : Composer les transformations

Multiplier des matrices revient à composer les transformations, c'est à dire les appliquer les unes après les autres (en partant de la plus intérieure): $B \circ C$ revient à appliquer la transformation $C$, puis $B$.

Exercice 4 :

Soient les matrices \begin{equation} C = \left [ \begin{array}{cc} 1 & 1 \\ 0 & 1 \\ \end{array} \right ] \text{ et } R = \left [ \begin{array}{cc} 0 & -1 \\ 1 & 0 \\ \end{array} \right ] \end{equation}

1) De quels types de transformations s'agit'il ?

2) L'opération matricielle est-elle commutative dans ce cas ?

# Correction exercice 4 :

C = np.array([[1,1],
              [0,1]]) # Cisaillement 
R = np.array([[0,-1],
             [1,0]]) # Rotation de 90 degrés.

fig, ax = plt.subplots(1,4,figsize=(24,6))
display_transform(C, ax[0], title='Cisaillement', name='C')
display_transform(R, ax[1], title='Rotation', name='R') 
display_transform(C@R, ax[2], title='Rotation puis cisaillement', name=r'C\circ R')
display_transform(R@C, ax[3], title='Cisaillement puis rotation', name=r'R\circ C');

plt.savefig('rotation.pdf')

Cela illustre pourquoi la multiplication de matrice n'est pas commutative en général. Appliquer la rotation puis le cisaillement n'est pas équivalent à appliquer le cisaillement puis la rotation. D'où :

$$BC \neq CB $$

Inverser une matrice : chercher la transformation inverse

Inverser une matrice revient à trouver la transformation qui annule la transformation associée à cette matrice. Si on garde en tête le fait que la multiplication est une composition des transformations (on en applique l'une puis l'autre) on comprend bien pourquoi: $$A^{-1}A = AA^{-1} = I$$

Exercice 5 :

1) Inverser les matrices A, C et R. Répresenter ces transformations inverses sur le carré unité. (Indice : utiliser np.linalg.inv).

2) Soit la matrice \begin{equation} D(t) = \left [ \begin{array}{cc} cos(t) & 0 \\ sin(t) & 1 \\ \end{array} \right ]\end{equation}

Regarder, pour plusieurs valeurs de t, les matrices D(t), leur inverse et leur déterminant, notamment quand le volume du carré unité passe par 0. Calculer la valeur du déterminant de chaque matrice. (Indice : utiliser np.linalg.det).

# Correction exercice 5 - question 1 :

invA = np.linalg.inv(A) # Inverse numériquement la matrice. 

fig, ax = plt.subplots(1,3,figsize=(24,8))
display_transform(A, ax[0])
display_transform(invA, ax[1], name='A^{-1}')
display_transform(A@invA, ax[2], name='AA^{-1}'); # Remarquez l'erreur numérique ce n'est pas exactement I

plt.savefig('inverse.pdf')


invC = np.linalg.inv(C) # Inverse numériquement la matrice. 

fig, ax = plt.subplots(1,3,figsize=(24,8))
display_transform(C, ax[0])
display_transform(invC, ax[1], name='C^{-1}')
display_transform(C@invC, ax[2], name='CC^{-1}')


invR = np.linalg.inv(R) # Inverse numériquement la matrice. 

fig, ax = plt.subplots(1,3,figsize=(24,8))
display_transform(R, ax[0])
display_transform(invR, ax[1], name='R^{-1}')
display_transform(R@invR, ax[2], name='RR^{-1}')

<matplotlib.axes._subplots.AxesSubplot at 0x11e943828>

# Correction exercice 5 - question 2 :


# Regardons quelques matrices et leur déterminant, quand le volume du carré unité passe par 0. 
fig, axes = plt.subplots(1,5, figsize=(5*6,6))
theta = np.linspace(0, np.pi)
j = 0
dots = []
detlist = []
for i,t in enumerate(theta): 
    # Cette famille de matrice ne change pas le second vecteur de la base canonique,
    # et fait décrire le cercle unité au premier.
    M = np.array([[np.cos(t),0],
                  [np.sin(t),1]])
    # np.linalg.det permet de calculer le déterminant. 
    detlist.append(np.linalg.det(M))
    
    if i>1 and i%(len(theta)//4) == 0:
        t = np.round(t, 2)
        M = np.array([[np.cos(t),0],
              [np.sin(t),1]])
        j+=1
        display_transform(M, axes[j], title='t={}'.format(t), name='M')
        dots.append((t,np.linalg.det(M)))
        
axes[0].plot(theta, detlist)
axes[0].scatter(*zip(*dots))
axes[0].hlines(0, theta.min(),theta.max())
axes[0].vlines(0,-1,1)

plt.savefig('non_inversible.pdf')

Cela illustre bien pourquoi il est impossible d'inverser une matrice dont le déterminant est nul, en effet, cela signifie que les volumes (ici les aires) sont "écrasés" en un ensemble de dimension inférieure (ici segment ou point). Il est impossible de trouver une transformation linéaire qui permette de retrouver le volume de départ à partir du segment.

Diagonaliser une matrice: Chercher la directions dans laquelle la transformation est une dilatation

Exercice 6 :

Diagonalisation :

1) Pour chacune des matrices suivantes, calculer les valeurs propres, répresenter la déformation du carré unité ainsi que les axes propres et vecteurs propres qui caractérisent ces transformations. (Indice : utiliser np.linalg.eig).

\begin{equation} I = \left [ \begin{array}{cc} 2 & 0 \\ 0 & 1 \\ \end{array} \right ] , J = \left [ \begin{array}{cc} 1 & 0 \\ 0 & 0.5 \\ \end{array} \right ] , K = \left [ \begin{array}{cc} 2 & 0 \\ 0 & 0 \\ \end{array} \right ] , L = \left [ \begin{array}{cc} 2 & 0 \\ 0 & -1 \\ \end{array} \right ] \text{ et } M = \left [ \begin{array}{cc} -1 & 0 \\ 0 & -2 \\ \end{array} \right ] \end{equation}\begin{equation} A = \left [ \begin{array}{cc} 1.2 & 0.46 \\ 0.45 & 0.8 \\ \end{array} \right ], B = \left [ \begin{array}{cc} 0.55 & 0.46 \\ 0.45 & 0.5 \\ \end{array} \right ], \text{ et } C = \left [ \begin{array}{cc} 0 & 0.9 \\ 1 & 0 \\ \end{array} \right ] \end{equation}

2) Superposer des trajectoires sur ces représentations, en proposant plusieurs conditions initiales

3) Interpréter, si possible, ces transformations en terme d'évolution du système Réticulum-Golgi.

# Correction exercice 6 - question 1 :


def display_eigenvectors(A, ax=None):
    """Display eigenvectors of the matrix A"""
    if ax is None:
        ax = plt.gca()

    # Calcule les vecteurs et valeurs propres ("eigenvectors", "eigenvalues")
    eigval, eigvect = np.linalg.eig(A)
    
    # Représentation les vecteurs propres
    # S'il s'agit de vecteurs propres complexes, ce n'est pas possible de les représenter de la sorte
    # et une erreur s'affiche.
    arrow ={'linewidth':2, 'head_width':0.05, 'length_includes_head':True}
    try:
        ax.arrow(0, 0, eigval[0]*eigvect[0,0], eigval[0]*eigvect[1,0], color='C1', label='e', **arrow)
        ax.arrow(0, 0,  eigval[1]*eigvect[0,1],eigval[1]* eigvect[1,1], color='C1', **arrow)
    except TypeError:
        pass
    
    # Représentation des axes propres
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    xmax = 10*np.max([ax.get_xlim(),ax.get_ylim()])
    try:
        ax.plot([0,eigvect[0,0]*xmax],[0,eigvect[1,0]*xmax], color='C1',alpha=.3)
        ax.plot([0,eigvect[0,1]*xmax],[0,eigvect[1,1]*xmax], color='C1',alpha=.3)
        ax.plot([0,-eigvect[0,0]*xmax],[0,-eigvect[1,0]*xmax], color='C1',alpha=.3)
        ax.plot([0,-eigvect[0,1]*xmax],[0,-eigvect[1,1]*xmax], color='C1',alpha=.3)
    except TypeError:
        pass
    
    ax.set(xlim=xlim, ylim=ylim)
    return ax


I = np.array([[2,0],
             [0,1]])
display_transform(I,name="I")
display_eigenvectors(I)


J = np.array([[1,0],
             [0,0.5]])
display_transform(J,name="J")
display_eigenvectors(J)


K = np.array([[2,0],
             [0,0.0001]])
display_transform(K,name="K")
display_eigenvectors(K)


L = np.array([[2,0],
             [0,-1]])
display_transform(L,name="L")
display_eigenvectors(L)


M = np.array([[-1,0],
             [0,-2]])
display_transform(M,name="M")
display_eigenvectors(M)


A = np.array([[1.2,0.46],
             [0.45,0.8]])
display_transform(A,name="A")
display_eigenvectors(A)
eigval, eigvect = np.linalg.eig(A)
eigval # 1.49699095, 0.50300905

B = np.array([[0.55,0.46],
             [0.45,0.5]])
display_transform(B,name="B")
display_eigenvectors(B)
eigval, eigvect = np.linalg.eig(B)
eigval # 0.98065886, 0.06934114


C = np.array([[0,0.9],
             [1,0]])
display_transform(C,name="C")
display_eigenvectors(C)
eigval, eigvect = np.linalg.eig(C)
eigval # 0.9486833, -0.9486833

array([ 0.9486833, -0.9486833])

# Correction exercice 6 - question 2 :

def display_trajectoire_ee(traj, ax):
    """Affiche la trajectoire `traj` dans l'espace d'état"""
    ax.plot(*zip(*traj), color='C3')
    ax.scatter(*zip(*traj), color='C3')


u0 = [0.1,.4]
traj = simulate(I,u0, T=100)
ax = display_transform(I,name="I")
display_eigenvectors(I)
display_trajectoire_ee(traj,ax)

u0 = [0.1,.4]
traj = simulate(J,u0, T=100)
ax = display_transform(J,name="J")
display_eigenvectors(J)
display_trajectoire_ee(traj,ax)

u0 = [0.1,.4]
traj = simulate(K,u0, T=100)
ax = display_transform(K,name="K")
display_eigenvectors(K)
display_trajectoire_ee(traj,ax)

u0 = [0.1,.4]
traj = simulate(M,u0, T=100)
ax = display_transform(M,name="M")
display_eigenvectors(M)
display_trajectoire_ee(traj,ax)

u0 = [0.1,.4]
traj = simulate(L,u0, T=100)
ax = display_transform(L,name="L")
display_eigenvectors(L)
display_trajectoire_ee(traj,ax)

u0 = [0.1,.4]
traj = simulate(A,u0, T=100)
ax = display_transform(A,name="A")
display_eigenvectors(A)
display_trajectoire_ee(traj,ax)

u0 = [0.1,.4]
traj = simulate(B,u0, T=100)
ax = display_transform(B,name="B")
display_eigenvectors(B)
display_trajectoire_ee(traj,ax)


u0 = [0.1,.4]
traj = simulate(C,u0, T=100)
ax = display_transform(C,name="C")
display_eigenvectors(C)
display_trajectoire_ee(traj,ax)

Perron-Frobenius

Exercice 7 :

Appliquer Perron-Frobenius

1) Pour chacune des matrices suivantes, calculer les valeurs propres, répresenter la déformation du carré unité ainsi que les axes propres et vecteurs propres qui caractérisent ces transformations.

\begin{equation} A = \left [ \begin{array}{cc} 1.2 & 0.46 \\ 0.45 & 0.8 \\ \end{array} \right ], B = \left [ \begin{array}{cc} 0.55 & 0.46 \\ 0.45 & 0.5 \\ \end{array} \right ], D = \left [ \begin{array}{cc} cos(\pi/4) & -sin(\pi/4) \\ sin(\pi/4) & cos(\pi/4) \\ \end{array} \right ], E = \left [ \begin{array}{cc} 0 & -1 \\ 1 & 0 \\ \end{array} \right ], F = \left [ \begin{array}{cc} 0 & 1 \\ 1.1 & 1.2 \\ \end{array} \right ], G = \left [ \begin{array}{cc} 1 & 1 \\ 0 & 1 \\ \end{array} \right ], H = \left [ \begin{array}{cc} 0 & 1 \\ 1 & 0 \\ \end{array} \right ], \text{ et } I = \left [ \begin{array}{cc} 0 & 1.05 \\ 1.06 & 0 \\ \end{array} \right ] \end{equation}

2) Quelles matrices sont positives ? Régulières ? Appliquer Perron-Frobenius.

3) Superposer des trajectoires sur ces représentations, en proposant plusieurs conditions initiales

4) Interpréter, si possible, ces transformations en terme d'évolution du système Réticulum-Golgi.

u0 = [0.1,.4]

A = np.array([[1.2,0.46],
             [0.45,0.8]])
traj = simulate(A,u0, T=100)
ax = display_transform(A,name="A")
display_eigenvectors(A)
display_trajectoire_ee(traj,ax)
eigval, eigvect = np.linalg.eig(A)
eigval # 1.49699095, 0.50300905


B = np.array([[0.55,0.46],
             [0.45,0.5]])
traj = simulate(B,u0, T=100)
ax = display_transform(B,name="B")
display_eigenvectors(B)
display_trajectoire_ee(traj,ax)
eigval, eigvect = np.linalg.eig(B)
eigval # 0.98065886, 0.06934114


# Rotation 
C= np.array([[np.cos(np.pi/4),-np.sin(np.pi/4)],
             [np.sin(np.pi/4),np.cos(np.pi/4)]])
traj = simulate(C,u0, T=100)
ax = display_transform(C,name="C")
display_eigenvectors(C)
display_trajectoire_ee(traj,ax)
eigval, eigvect = np.linalg.eig(C)
eigval # 0.70710678+0.70710678j, 0.70710678-0.70710678j


# Rotation                    
D= np.array([[0,0.95],
             [-0.95,0]])
traj = simulate(D,u0, T=100)
ax = display_transform(D,name="D")
display_eigenvectors(D)
display_trajectoire_ee(traj,ax)
eigval, eigvect = np.linalg.eig(D)
eigval # 0.+0.95j, 0.-0.95j

# Rotation                    
E= np.array([[0,+1.05],
             [-1.05,0]])
traj = simulate(E,u0, T=100)
ax = display_transform(E,name="E")
display_eigenvectors(E)
display_trajectoire_ee(traj,ax)
eigval, eigvect = np.linalg.eig(E)
eigval # 0.+1.05j, 0.-1.05j

# Les matrices de rotation ne vérifient pas l'hypothèse de positivité de Perron-Frobenius.


# Regular
F= np.array([[0,1],
             [1.1,1.2]])
traj = simulate(F,u0, T=100)
ax = display_transform(F,name="F")
display_eigenvectors(F)
display_trajectoire_ee(traj,ax)
eigval, eigvect = np.linalg.eig(F)
eigval # -0.6083046,  1.8083046       
    
# Cisaillement
G=np.array([[1,1],
            [0,1]])
traj = simulate(G,u0, T=100)
ax = display_transform(G,name="G")
display_eigenvectors(G)
display_trajectoire_ee(traj,ax)
                     
# Flip
H = np.array([[0,1],
              [1,0]])
traj = simulate(H,u0, T=100)
ax = display_transform(H,name="H")
display_eigenvectors(H)
display_trajectoire_ee(traj,ax)
eigval, eigvect = np.linalg.eig(H)
eigval # 1., -1.

# Oscillante
I =np.array([[0,1.05],
             [1.06,0]])
traj = simulate(I,u0, T=100)
ax = display_transform(I,name="I")
display_eigenvectors(I)
display_trajectoire_ee(traj,ax)
eigval, eigvect = np.linalg.eig(I)
eigval # 1.05498815, -1.05498815

# La matrices flip et oscillante ne sont pas régulières

/Users/bperez/miniconda3/lib/python3.7/site-packages/numpy/core/numeric.py:538: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

array([ 1.05498815, -1.05498815])

Partie 3 : Inférence

Exercice 8 :

Une protéine marquée à la GFP est suivie au cours du temps dans les cellules pancréatiques et les concentrations suivantes ont été mesurées dans le Reticulum et le Golgi.

temps	départ	5 minutes	10 minutes	15 minutes	20 minutes
Reticulum	1.0	1.2	1.3	1.4	1.6
Golgi	0.0	0.3	0.5	0.7	0.8

1) Modéliser ce flux protéique à l'aide d'un modèle linéaire. En utilisant la méthode des moindres carrés, estimer les valeurs des paramètres du modèle linéaire décrivant le mieux le comportement du système.

2) Quelles concentrations prédisez-vous 40 minutes après le début de l'expérience ? Après 5 heures ?

# Correction Exercice 8 - Question 1 

# Données 

data=np.array([[1.0,1.2,1.3,1.5,1.5],
               [0.0,0.2,0.5,0.7,0.8]])

# Modèle linéaire et déviation par rapport aux données

def erreur_modele(a,b,c,d,data):
    A=np.array([[a,b],
                [c,d]])
    
    traj=[[1,0],]
    erreur=0
    for i in range(4):
        traj.append(A@traj[-1])
        erreur=erreur+sum((data[:,i+1] - traj[-1])**2)
    return(erreur)

a=1.1
b=0.1
c=0.3
d=0.6

erreur=erreur_modele(a,b,c,d,data)

# Inférence des paramètres - par méthode des moindres carrés 
minium_erreur=10000
best_a=10
best_b=10
best_c=10
best_d=10
for a in np.linspace(0,2,21):
    for b in np.linspace(0,2,21):
        for c in np.linspace(0,2,21):
            for d in np.linspace(0,2,21):
                erreur=erreur_modele(a,b,c,d,data)
                if erreur<minium_erreur:
                    minium_erreur=erreur
                    best_a=a
                    best_b=b
                    best_c=c
                    best_d=d


print("a: {}, b: {}, c: {}, d: {}, erreur {}".format(round(best_a,1),round(best_b,1),round(best_c,1),round(best_d,1), minium_erreur))

a: 1.1, b: 0.1, c: 0.3, d: 0.6, erreur 0.04783918000000006