Séance 3 : Simuler l'ouverture de canaux ioniques grâce aux chaînes de Markov¶
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
%matplotlib inline
Modélisation¶
Nous allons utiliser une Chaîne de Markov pour modéliser l'ouverture d'un canal ionique :
- avec comme espace d'état (de cardinal fini) l'ensemble
{Fermé, Ouvert}ou $\{0,1\}$. - avec comme condition initiale $X_0$, qui peut être une variable aléatoire (Sa loi est appelée mesure initiale de la chaîne) ou une valeur déterministe (par exemple
Ouvert, ce qui revient à avoir comme mesure initiale $P(X_0=1)=1$). - avec les probabilités de transitions: $a_{ij} = \mathbb P (X_{n+1} = j \mid X_n=i)$
On note la matrice de transition suivante: $$ A = \begin{pmatrix} 1-p & p \\ q & 1-q \\ \end{pmatrix} $$
Exercice 1 :¶
1) Quelles sont les propriétés d'une matrice de transition si $p \notin \{0,1\} $ and $q \notin \{0,1\}$ ?
2) Écrire une fonction qui prend une matrice comme paramètre et lève une erreur si ce n'est pas une matrice stochastique. (Indice: Utiliser assert, numpy.testing.assert_allclose).
## Réaliser des tirages de nombres aléatoires à l'aide d'un ordinateur.
# Fixe une graine
np.random.seed(111)
# Génère un nombre pseudo aléatoire
assert np.random.random() - 0.6121701756176187 < 1e-10
Exercice 2 :¶
1) Écrire une fonction mc_trajectoire(mu, P, Tmax) qui échantillonne, du temps $0$ à $Tmax$, la chaîne de Markov ($mu$, $P$) où $P$ est une matrice de transition 2x2 et $mu$ la mesure initiale. (Indice : Utiliser np.random.random).
2) Représenter la trajectoire de la chaîne (fonction trace_mc_trajectoire) et appliquer aux chaînes suivantes :
P = np.array([[0.9,.1],
[0.4,0.6]])
P3 = np.array([[0,1],
[1,0]])
mu = np.array([.3,.7])
mu3 = np.array([.1,.9])
B) Rajout d'un troisième état :¶
Exercice 3 :¶
1) Écrire une fonction sample_cdf qui permet d'échantillonner une loi de probabilité discrète arbitraire. (Indice : Utiliser np.random, np.cumsum, np.sum).
2) Ré-écrire la fonction mc_trajectoire pour prendre en compte des matrices stochastique de taille arbitraire.
Utiliser votre fonction sample_cdf ou np.random.choice.
L'appliquer pour la chaîne suivante :
P2 = np.array([[0.8,.1,.1],
[0.4,0.3,0.3],
[0.2,0.6,0.2]])
mu2 = np.array([.3,.7,0])
Partie II : Mesure invariante et théorème ergodique :¶
Exercice 4 :¶
Écrire une fonction qui renvoie la mesure invariante d'une matrice stochastique $P$. (Indice : Utiliser np.linalg.eig et np.argmax).
A) Proportion du temps passé en un état donné :¶
Exercice 5 :¶
1) Écrire une fonction trace_proportion_temps_passe qui affiche la proportion du temps passé dans chaque état pour une trajectoire.
L'appliquer sur la matrice $P2$. Tester avec 100 trajectoires de 50 pas de temps, puis 1000 trajectoires de 500 pas de temps.
Vérifier la convergence vers la mesure invariante.
(Indice : Utiliser plt.hist et plt.hlines).
2) Représenter la convergence vers la mesure invariante au cours du temps.
B) Probabilité d'être dans un état donné à un temps donné :¶
Exercice 6 :¶
1) Écrire une fonction qui estime numériquement $P(X_t)$ à partir de plusieurs trajectoires de $P2$.
2) Que se passe-t-il pour $P3$ ? Comment l'expliquer ?
III) Inférence :¶
Exercice 7 :¶
Les fichiers ion_channel_1 et ion_channel_2 contiennent des séries temporelles de mesure de patch clamp qui donnent l'état ouvert ou fermé d'un canal au cours du temps.
Écrire une fonction estimate_P(trajectoire) qui estime une matrice de transition de taille 2x2 à partir de la trajectoire.
(Indice : Utiliser np.fromfile et zip).
Exercice 8 :¶
1) Écrire une fonction qui renvoie les liste des durées pendant lequel le canal reste ouvert (resp. fermé).
Tracer les deux histogrames pour ion_channel_1 (ouvert/fermé) et ion_channel_2 (ouvert/fermé). (Indice : Utiliser plt.hist).
2) Quelle est la distribution théorique de ces temps d'ouverture/fermeture si le phénomène suit la distribution d'une chaîne de Markov de taille 2x2 ? La tracer par dessus les histogrammes.
(Indice : Utiliser: scipy.stats.xxx.pmf où xxx est le nom de la loi discrète parmi: https://docs.scipy.org/doc/scipy/reference/stats.html#discrete-distributions )
Les deux jeux de données donnent la même matrice de transition mais le second n'est pas issu d'une chaîne de Markov à deux états ! La méthode pour étudier ce genre de sytème est de faire un modèle de mélange (statistique) pour voir si le temps de séjour n'est pas la composition de plusieurs loi géométriques, ce qui donne un indice sur le nombre d'états cachés qui existent. On peut ensuite paramétriser un modèle de Markov Caché.
# Code pour générer le jeu de données...
P2 = np.array([[0.8,.1,.1],
[0.9,0,0.1],
[0.2,0,.8]])
mu2 = np.array([.3,.7,0])
def assert_transition_matrix(P):
"""Raise an Assertion error if the matrix is not stochastic"""
assert P.ndim == 2, 'Not a 2D matrix'
assert P.shape[0] == P.shape[1], 'Matrix is not square'
np.testing.assert_allclose(np.sum(P,1), [1]*P.shape[0], err_msg='Lines do not sum at 1.')
def mc_trajectoire(mu, P, Tmax):
'''Simule une chaîne de markov discrète de matrice de transition P (de dimension finie)
et de mesure initiale mu sur Tmax pas de temps'''
# Tester si la matrice de transition est valide.
assert_transition_matrix(P)
assert mu.sum() - 1 < 1e-6, 'mesure initiale invalide'
# Initialisation
traj = np.zeros(Tmax, dtype=int)
choices = np.arange(P.shape[0], dtype=int)
traj[0] = np.random.choice(choices, 1, p=mu)
# Simulation
for i in range(Tmax-1):
traj[i+1] = np.random.choice(choices,1,p=P[traj[i],:])
return traj
def trouver_mesure_invariante(P):
"""Retourne la mesure invariante d'une matrice stochastique P"""
# Tester si la matrice est stochastique.
assert_transition_matrix(P)
# Trouver le vecteur propre à gauche
# /!\ linalg.eig renvoie la valeur propre à droite,
# il faut transposer la matrice de départ pour avoir
# celle à gauche.
valp, vecp = np.linalg.eig(P.transpose())
# Trouver la valeur propre dominante
# On sait qu'elle existe graçe au théorème de Perron-Frobenius.
pos_valp_dom = np.argmax(valp)
# Trouver le vecteur propre dominant
vecp_dom = vecp[:, pos_valp_dom]
# Vérifier que la valeur propre dominante est 1.
assert abs(np.max(valp) - 1) < 1e-10
# Remettre le vecteur propre dominant à l'échelle pour que la somme
# de ses éléments soit 1.
vecp_dom = vecp_dom/np.sum(vecp_dom)
return vecp_dom
def collapse(P):
o_to_c = P[0,1] + P[0,2]
pi = trouver_mesure_invariante(P)
c_to_o = (pi[1] * P[1,0] + pi[2] * P[2,0]) / (pi[1]+pi[2])
return np.array([[1-o_to_c, o_to_c],
[c_to_o, 1-c_to_o]])
P1 = collapse(P2)
print(P1)
ion_channel_2 = mc_trajectoire(mu2, P2, 100000)
ion_channel_1 = mc_trajectoire(trouver_mesure_invariante(P1), P1, 100000)
ion_channel_2 = np.int_(np.bool_(ion_channel_2))
ion_channel_1.tofile('ion_channel_1', sep=',')
ion_channel_2.tofile('ion_channel_2', sep=',')