Chaînes de Markov et Canaux Ioniques¶
# Commençons par importer les outils dont nous aurons besoin.
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
%matplotlib inline
Une histoire de canal¶
Soit un transporteur transmembranaire que l'on trouve dans deux conformations: ouvert ou fermé. On sait que la protéine change de conformation de manière aléatoire selon un certain taux que l'on peut mesurer.
On souhaite modéliser ce système pour répondre à un certain nombre de questions:
- Quel est le comportement d'un tel transporteur sur le temps long ?
- Quel est le comportement d'une population de transporteurs ?
En utilisant les outils du cours, quelle approche utiliseriez vous pour modéliser ce système?
Exercice: Écrivez un tel modèle
Modélisation¶
Nous alons utiliser une Chaîne de Markov
- 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)$
Que l'on note la matrice de transition suivante: $$ A = \begin{pmatrix} 1-p & p \\ q & 1-q \\ \end{pmatrix} $$
Si $p \notin \{0,1\} $ and $q \notin \{0,1\}$, Cette chaîne de markov est irréductible car elle comporte une seule classe d'état qui communiquent. Cette classe est récurrente. Comme la chaîne a un espace d'état fini, cette récurrence est nécessairement positive.
Quelles sont les propriétés d'une matrice de transition ?
Exercice: Écrivez une fonction en python qui prend une matrice comme paramètre et lève une erreur si ce n'est pas une matrice stochastique.
P = np.array([[0.9,.1],
[0.4,0.6]])
P2 = np.array([[0.8,.1,.1],
[0.4,0.3,0.3],
[0.2,0.6,0.2]])
P3 = np.array([[0,1],
[1,0]])
mu = np.array([.3,.7])
mu3 = np.array([.1,.9])
mu2 = np.array([.3,.7,0])
# Écrire une fonction assert_transition_matrix(P)
# qui lève une exception si la matrice P
# N'est pas une matrice stochastique.
# Utilisez: assert, np.testing.assert_allclose
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.')
P = np.array([[0.9,.1],[.2,.8]])
assert_transition_matrix(P)
# Essayez:
#P = np.array([[0.8,.1],[.2,.8]])
#assert_transition_matrix(P)
Remarque: Comparaison de nombre à virgule flotante¶
Vous observez que l'on ne teste (jamais) la somme de nombres à virgule flotante avec une égalité stricte. En effet, la représentation des nombres à virgule flottante pose de nombreux problèmes.
Essayez:
assert 0.1 + 0.2 == 0.3
Vous devriez avoir une erreur. En effet s'il est impossible de représenter $\frac{1}{3}$ dans l'écriture décimale sans utiliser une infinité de nombres après la virgule ($0.33333...$), c'est aussi le cas de $0.1$ en écriture binaire ($0.0001100110011...$)
| Nombre | Base 10 | Base 2 |
|---|---|---|
| 1 | 1 | 1 |
| 2 | 2 | 10 |
| 3 | 3 | 11 |
| $\frac{1}{2}$ | 0.5 | 0.1 |
| $\frac{1}{10}$ | 0.1 | $0.0\overline{0011}$ |
| $\frac{1}{5}$ | 0.2 | $0.\overline{0011}$ |
| $\frac{1}{3}$ | $0.\overline{3}$ | $0.\overline{01}$ |
Il y a 10 types de personnes dans le monde: ceux qui comprennent le binaire, ceux qui ne le comprennent pas et ceux qui pensaient que cette blague était en ternaire.
En général on teste l'égalité des nombre à virgule flotante à une certaine précision près. Par exemple on testera:
assert abs(((0.1 + 0.2)-0.3)/0.3) < 1e-6
La fonction numpy.testing.assert_allclose permet de faire de telles comparaisons élément par élément sur un array.
Une seconde solution est d'utiliser une bibliothèque de calcul avec des nombre décimaux à précision arbitraire:
from decimal import Decimal
assert Decimal('0.2')+Decimal('0.1') == Decimal('0.3')
Mais ces bibliothèques sont beaucoup plus lentes que les calculs natifs que fait le processeur en virgule flotante. Pour la majeure partie des applications en biologie, cela n'est pas nécessaire.
References :
- Floating point guide
- Goldberg, David (March 1991). "What Every Computer Scientist Should Know About Floating-Point Arithmetic" ACM Computing Surveys. 23 (1): 5–48. doi:10.1145/103162.103163
Simulation¶
Comment simuler la trajectoire d'un tel objet ?
La trajectoire est la donnée $(X_n)_{n = 0,1,2...}$, c'est à dire une suite à valeur dans l'espace d'état et indexée sur les nombre entiers.
Réaliser des tirages de nombres aléatoires à l'aide d'un ordinateur.¶
En fait, il est impossible de réaliser des tirages de nombre aléatoire avec la majorité des ordinateurs sur lequels vous serez amené à travailler. Il faudrait un matériel particulier. À la place on utilise des suites de nombres pseudo aléatoires.
À la place les ordinateurs utilisent un algorithme qui génère une suite de nombres, qui a l'air aléatoire. Par exemple un générateur congruentiel linéaire fonctionne ainsi:
\begin{equation} X_{n+1} = (a X_n + b) \text{ mod } m \end{equation}
Cette suite est périodique et si on prend a = 25, c = 16, m = 256, X0 = 50 cette periode est de 5 !
50, 242, 178, 114, 50, 242, 178, 114, 50, 242, ...
Il est donc important de bien choisir ces paramètres.
Numpy implémente un générateur de nombre pseudo aléatoire appelé Mersenne Twister, qui est très rapide et génère une chaine de nombres en 32bits uniforméments distribués. Chaque appel à la fonction np.random.random() va vous donner le prochain élément dans la suite de nombres pseudo aléatoires qu'utilise numpy. La période de cette suite est de $2^{{19937}}-1$.
Les nombres pseudo aléatoires dépendent de la condition initiale du générateur- que l'on appelle la graine. Vous pouvez fixer la graine de numpy à l'aide de la commande np.seed. Essayez donc:
np.random.seed(111)
# Pas d'erreur ! On peut lancer cette cellule autant de fois que l'on veut !
assert np.random.random() - 0.6121701756176187 < 1e-10
Méthode naïve¶
$a_{ij} = \mathbb P (X_{n+1} = j \mid X_n=i)$
# Écrire une fonction mc_trajectoire(mu, P, Tmax) qui
# échantillonne la chaîne de markov (mu, P) où P est une matrice 2x2 du temps 0 au temps Tmax.
# Utilisez: np.random
def simulation_naive():
pass
def mc_trajctoire(mu, P, Tmax):
"""Retourne la trajectoire de la chaine mu,P sur 1...Tmax"""
trajectoire = []
# Condition initiales:
x0 = int(np.random.random() < mu[1])
trajectoires.append(x0)
# Chaque pas de temps
for t in range(Tmax-1):
x = trajectoire[-1]
rd = np.random.random()
if x == 0:
if rd < P[0,1]:
new_x = 1
else:
new_x = 0
else:
new_x = int(rd < P[1,1])
trajectoire.append(new_x)
return trajectoire
Si on prend une chaine légèrement plus compliquée, par exemple en ajoutant un troisième état, il faut changer notre algorithme !
Peut-on trouver un algorithme général ?
# (*) Écrire une fonction choice qui permet d'échantillonner une
# loi de probabilité discrète arbitraire
# Utilisez: np.random, np.cumsum, np.sum
Comment Échantillonner une densité de probabilité discrète ?¶
Le problème principal de cet algorithme est de choisir l'état au temps $n+1$ sachant que la chaine était dans l'état $i$ au temps $n$. Cela revient à échantillonner la densité de probabilité discrète définie sur la $i$-ème ligne de la matrice de transition.
Le problème c'est qu'un ordinateur par défaut ne peux générer que des nombres pseudo-aléatoires entre 0 et 1 (np.random.random()). Comment faire ?
Méthode de la fonction de répartion inverse¶
Cette méthode fait appel à l'inverse de fonction de répartition. La fonction de répartition $F$ donne $F(k) = P(X\leq k)$. Elle est constante par morceaux pour une variable discrète. Il suffit de l'inverser, de tirer un nombre entre 0 et 1 et de regarder sur quel "palier" on tombe pour avoir la valeur recherchée.
def sample_cdf(p):
''' Échantillonne une distribution discrète.
Args:
p (np.array de longueur 4): p[i] probabilités de X=i (pas necessairement normalisées pour sommer à 1)
Return:
i \in 0...3. '''
# Normaliser si necessaire.
if p.sum() =! 1:
p /= p.sum()
# Tirage.
rdm = np.random.random()
for i,csum in enumerate(np.cumsum(p)):
if rdm < csum:
return i
En pratique¶
Pour plus d'efficacité vous pouvez utiliser la fonction np.random.choice(<espace d'etat>, <nombre de tirage>, p=<probabilités>). Comme numpy est libre et Open-source, vous pouvez aller vérifier que cette fonction utilise en fait la méthode de la fonction de répartition en triant les valeurs pour minimiser le nombre d'opérations.
Aller plus loin...¶
Exercice: Réfléchir à l'efficacité de cette méthode s'il y a un grand nombre d'états possibles ou si les états ont des probabilités très différentes. Quand est-il de la méthode naïve par dichotomie ?
Pour la réponse et probablement plus de détails que vous n'en aurez jamais besoin, les passionnés pourront se tourner vers Darts, Dice, and Coins: Sampling from a Discrete Distribution, une très intéressante comparaison des algorithmes avec leurs coûts en temps d'initialisation, de calcul et en mémoire dans les meilleurs et pire des cas.
# Ré-écrivez la fonction mc_trajectoire pour prendre en compte des matrices stochastique de taille arbitraire.
# Utilisez votre fonction choice ou np.random.choice.
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 trace_mc_trajectoire(traj, noms_etats=None, ax=None):
"""Trace la trajectoire d'une chaîne de markov.
Args:
traj (np.array): traj[i] est l'état de la chaîne à l'instant i.
noms_etats (iterable): nom associé à chaque état
ax (matplotlib.axis): axis on which the trajectory is drawn
"""
if ax is None:
fig, ax = plt.subplots(figsize=(15,2))
ticks = np.arange(traj.max()+1)
if noms_etats is None:
noms_etats = ticks
ax.plot(np.arange(len(traj)), traj,color='k')
ax.set(xlabel='Temps', ylabel='État', yticks=ticks, yticklabels=noms_etats)
ax.scatter(np.arange(len(traj)), traj)
P = np.array([[0.9,.1],
[0.4,0.6]])
mu = np.array([.3,.7])
traj = mc_trajectoire(mu, P, 100)
trace_mc_trajectoire(traj, ['Fermé','Ouvert'])
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])
traj = mc_trajectoire(mu2, P2, 100)
trace_mc_trajectoire(traj, ['Fermé','Ouvert','Inactivé'])
P3 = np.array([[0,1],
[1,0]])
mu3 = np.array([.1,.9])
traj = mc_trajectoire(mu3, P3, 100)
trace_mc_trajectoire(traj, ['Fermé','Ouvert','Inactivé'])
Mesure invariante et théorème ergodique¶
La mesure invariante d'une chaine de Markov recurrente positive est donnée par le vecteur propre dominant de sa matrice de transition.
Si on a une trajectoire, la proportion du temps passé en $i$ est:
$$\tau_i = \frac{1}{T}\sum^T_{k=0} 1_{(X_k = i)}$$
Si on a N trajectoires notées $X^{(1)},X^{(2)},...X^{(N)}$, l'estimateur de la probabilité d'un evenement:
$$ \mathbb P (X_t = i) \approx \rho_{t,i} = \frac{1}{N} \sum_{j=0}^N 1_{(X^{(j)}_t = i)} $$
# Écrire une fonction trouver_mesure_invariante(P)
# Qui renvoie la mesure invariante d'une matrice stochastique P
# Utilisez: np.linalg.eig np.argmax
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. Exercice: le prouver.
valp, vecp = np.linalg.eig(P.transpose())
# Trouver la valeur propre dominante
# On sait qu'elle existe grâce au thèorème de Perron-Frobenius.
pos_valp_dom = np.argmax(valp)
vecp_dom = vecp[:, pos_valp_dom]
# Vérifier que la valeur propre dominante est 1.
assert abs(np.max(valp) - 1) < 1e-10
# Trouver le vecteur propre dominant
# et le remettre à l'échelle pour que la somme
# de ses éléments soit 1.
vecp_dom = vecp[:,pos_valp_dom]/np.sum(vecp_dom)
return vecp_dom
mesure_inv = trouver_mesure_invariante(P2)
np.testing.assert_allclose(mesure_inv, mesure_inv @ P2)
# Faire de nombreuses réalisations,
R = 10000 #Nombre de trajectoires
T = 500 #Longueur d'une trajectoire
state_space = np.arange(P2.shape[0])
names = ('Fermé', 'Ouvert', 'Inactivé')
mu = np.random.random(P2.shape[0])
mu /= mu.sum()
trajs = np.array([mc_trajectoire(mu,P2,T) for _ in range(R)])
print(trajs.shape)
mesure_inv = trouver_mesure_invariante(P2)
mesure_inv
Proportion du temps passé en un état donné¶
On cherche à mesurer la proportion du temps passé dans un état donné et à l'afficher sous forme graphique
# Écrire une fonction proportion_temps_passe qui affiche la proportion
# du temps passé dans chaque état pour une trajectoire.
# Vérifier la convergence vers la mesure invariante.
# Utilisez plt.hist, plt.hlines
def trace_proportion_temps_passe(espace_detat, trajs, mesure_inv=None, noms_etats=None, ax=None):
if ax is None:
fig, ax = plt.subplots(figsize=(7,5))
ticks = espace_detat
if noms_etats is None:
noms_etats = espace_detat
plt.hist(trajs.flat, density=True, bins=[0,1,2,3], label = r'$\frac{1}{T}\sum^T 1_{(X_k = i)}$')
if mesure_inv is not None:
ax.hlines(mesure_inv,
0,len(state_space),
label='Mesure invariante')
ax.set(xlabel='État',
ylabel='Proportion du temps passé',
title=r'$\frac{1}{T}\sum^T 1_{(X_k = i)}$',
xticks=state_space,
xticklabels=noms_etats)
plt.legend()
proportions = np.array([ np.sum(trajs==i) for i in state_space]) /trajs.size
trace_proportion_temps_passe(state_space,
trajs,
mesure_inv=mesure_inv,
noms_etats=names)
def trace_convergence_temps_passe(trajs, noms, mesure_inv, mu,ax=None):
proportions = np.array([np.array([(trajs[0,:t]==i).sum()/(trajs[0,:t]).size for i in state_space])
for t in range(1,trajs.shape[1])])
if ax is None:
fig, ax = plt.subplots(figsize=(16,4))
for i,name in enumerate(names):
plt.plot(proportions[:,i],label=name)
ax.hlines(mesure_inv,
0,trajs.shape[1],
linestyle=':',
color=('C0','C1','C2'),
label='Mesure invariante')
ax.hlines(mu,
0,trajs.shape[1],
linestyle='-.',
linewidth=1,
color=('C0','C1','C2'),
label='Mesure initiale')
ax.set(xlabel='Temps',
ylabel='Proportion du temps passé',
ylim=(0,1),
xlim=(0,trajs.shape[1]),
title=r'Convergence de $\frac{1}{T}\sum^T 1_{(X_k = i)}$')
ax.legend()
trace_convergence_temps_passe(trajs,names, mesure_inv, mu)
Probabilité d'être dans un état donné à un temps donné¶
# (*) Écrire une fonction qui estime numériquement P(X_t) à partir de plusieurs trajectoires.
# (**) Écrire une fonction qui calcule P(X_t) par récurrence.
# Vérifier la convergence en loi
def trace_convergence_loi(trajs, noms, mesure_inv, mu, ax=None):
proportions = np.array([(trajs==n).sum(0)/trajs.shape[0] for n,name in enumerate(noms)])
if ax is None:
fig, ax = plt.subplots(figsize=(16,4))
for i,name in enumerate(names):
plt.plot(proportions[i,:],label=name)
ax.hlines(mesure_inv,
0,trajs.shape[1],
linestyle=':',
color=('C0','C1','C2'),
label='Mesure invariante')
ax.hlines(mu,
0,trajs.shape[1],
linestyle='-.',
linewidth=1,
color=('C0','C1','C2'),
label='Mesure initiale')
ax.set(xlabel='Temps',
ylabel='Probabilité',
ylim=(0,1),
xlim=(0,trajs.shape[1]),
title=r'Convergence de $P(X_n = i)$')
ax.legend()
trace_convergence_loi(trajs[:100,:],names, mesure_inv, mu)
# Que se passe-t-il pour P3 ? Comment l'expliquer ?
# Faire de nombreuses réalisations,
R = 1000 # Nombre de trajectoires
T = 100 #Longueur d'une trajectoire
state_space = np.arange(P3.shape[0])
names = ('Fermé', 'Ouvert')
trajs = np.array([mc_trajectoire(mu3,P3,T) for _ in range(R)])
mesure_inv3 = trouver_mesure_invariante(P3)
trace_convergence_temps_passe(trajs,names,mesure_inv3,mu3)
plt.show()
trace_convergence_loi(trajs,names,mesure_inv3,mu3)
plt.show()
Retour à nos canaux¶
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
# Que constatez vous ?
def estimate_P(traj):
"""Estime la matrice de transition de la chaine 2x2 à partir d'une réalisation de la trajectoire.
Args:
traj (np.array): réalisation de la trajectoire.
Returns:
(np.array) Une matrice stochastique 2x2 .
"""
def estimate_P(traj):
"""Estime la matrice de transition de la chaine 2x2 à partir d'une réalisation de la trajectoire.
Args:
traj (np.array): réalisation de la trajectoire.
Returns:
(np.array) Une matrice stochastique 2x2 .
"""
# On vérifie que les données sont conformes.
assert all(np.logical_or(traj==0, traj==1))
# Créer une matrice vide.
P = np.zeros((2,2))
# Compter le nombre de transition entre chaque
# paire d'état
for pn,pnm1 in zip(traj[1:], traj[:-1]):
P[pnm1,pn] += 1
# Normaliser pour que la matrice soit stochastique.
P /= np.hstack([P.sum(1)[:,np.newaxis]]*2)
return P
ion_channel_1 = np.fromfile('ion_channel_1', sep=',', dtype=int)
ion_channel_2 = np.fromfile('ion_channel_2', sep=',', dtype=int)
fig, ax = plt.subplots(2,1,figsize=(20,3))
trace_mc_trajectoire(ion_channel_1[:100], ['Fermé','Ouvert'], ax[0])
trace_mc_trajectoire(ion_channel_2[:100], ['Fermé','Ouvert'], ax[1])
estimate_P(ion_channel_1)
estimate_P(ion_channel_2)
Comment savoir si notre modèle est juste ?
Dans un modèle à deux états, les temps de séjour dans l'état ouvert ouf fermé suivent lois géométriques. Regardons si c'est le cas.
# Écrire une fonction temps_changement(trajectoire) 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é).
# Utilisez: plt.hist,
# Quelle est la distribution théorique de ces temps d'ouverture/fermeture
# si le phénomène suit la distribution d'une chaine de markov de taille 2x2 ?
# Tracez-là par dessus les histogrammes.
# Utilisez: 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
def measure_residence_times(traj):
"""Estime les temps de séjour de la chaine 2x2 à partir d'une réalisation de la trajectoire.
Args:
traj (np.array): réalisation de la trajectoire.
Returns:
(tuple of list): liste des temps de séjour dans l'état 0, liste des temps de séjour dans l'état 1.
"""
# Vérifie qu'il s'agit bien d'une trajectoire de chaine à deux état.
assert all(np.logical_or(traj==0, traj==1)), 'Invalid trajectory'
# Initialisation
stretch = ([], [])
current = 1
# Pour chaque transition:
# si l'état a changé on enregistre la durée du séjour dans l'état précédent,
# si non on ajoute 1 à la durée du séjour.
for pn,pnm1 in zip(traj[1:], traj[:-1]):
if pn != pnm1:
stretch[pnm1].append(current)
current = 1
else:
current += 1
stretch[traj[-1]].append(current)
return stretch
assert measure_residence_times(np.array([0,0,0,1,1,1,0,0,1,0,0])) == ([3,2,2],[3,1])
def display_residence_times(rtimes):
"""Affiche le résultat de measure_residence_time"""
fig, axes = plt.subplots(1,2, figsize=(12,4))
for i,(val,ax) in enumerate(zip(rtimes,axes)):
x = np.arange(np.max(val))
ax.step(x+1, scipy.stats.geom.pmf(x, 1/np.mean(val)),label='Loi géométrique')
h,b,_ = ax.hist(val, density=True, bins=x, label='Empirique')
ax.set(xlabel='durée', ylabel='densité', title="Temps de séjour dans l'état {}".format(i))
ax.legend()
display_residence_times(measure_residence_times(ion_channel_1))
display_residence_times(measure_residence_times(ion_channel_2))
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 état ! Dans la littérature, la méthode pour étudier ce genre de sytème et 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ée...
P2 = np.array([[0.8,.1,.1],
[0.9,0,0.1],
[0.2,0,.8]])
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)
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=',')