import matplotlib.pyplot as plt
import numpy as np
import pickle
from scipy.optimize import curve_fit
w=np.pi/4 #Periode premiere phase
x=np.arange(0,20,0.0005) #Position x
def ratchet(x,maxx,w):
"""
:param x Array position:
:param maxx Maximum de la fonction:
:param w Periode phase ascendente du potentiel:
:return Potentiel en dent de scie:
"""
#W=w*((b+a)/a) #Periode seconde phase
W=1
a=maxx/w
b=-maxx/(W-w)
if x%W<=w:
return a*(x%W)
else:
return b*((x%W)-w)+a*w
print(w+w*((1+0.2)/0.2))
l=[]
for i in range(0,len(x)): #Boucle associe a chaque point de x une valeur de potentiel dans une liste
l.append(ratchet(x[i],1,1/3))
l=np.array(l) #Transforme la liste en array
t=np.arange(0,20,0.0005) #Array temps
def gauss(x):
"""
:param x Nombre d'itérations:
:return Une distribution normale de forces aléatoires:
"""
l=[]
for i in range(0,x):
l.append(np.random.normal(loc=0,scale=1))
return l
def vprime(x,maxx,w):
"""
:param x Array position:
:param maxx Maximum de la fonction:
:param w Periode phase ascendente du potentiel:
:return Dérivée du potentiel en dent de scie:
"""
W=1
a=maxx/w
b=-maxx/(W-w)
if x%W<=w:
return a
else:
return b
def periode(w):
"""
:param w Période
:return Array de 0 et 1 périodique
"""
x= np.sin((2*np.pi/w)*t)
y=np.zeros(len(x))
y[x<0]=0
y[x>=0]=1
return y
def euler(w,amp):
"""
:param w Période
:param amp Amplitude du potentiel
:return Solution de l'équation différentielle stochastique
"""
x=0
l=[]
for i in range(0,len(t)):
l.append(x)
x+=(np.random.normal(0,4.5)-vprime(x,amp,1/3)*periode(w)[i])*step
return l,l[-1]
t=np.arange(0,100,0.01) #Temps
step=0.01 #Pas temporel
x=np.arange(0,20,0.0005) #Array position x
vprim=[] #Liste vide dérivée potentiel
v=[] #Liste vite potentiel
for i in range(0,10000): #Ajoute chaques points du potentiel et de sa dérivée dans leurs listes respectives
vprim.append(vprime(x[i],10,1/3))
v.append(ratchet(x[i],10,1/3))
##PLOT
fig,ax=plt.subplots(figsize=(13,7))
plt.grid()
plt.plot(x[:10000],vprim,'--r',label='grad(Potentiel)') #Plot dérivée potentiel asymétrique
plt.plot(x[:10000],v,'k',label='Potentiel') #Plot potentiel asymétrique
ax.set_xlabel('Position sur x')
ax.set_ylabel('Potentiel')
plt.legend()
plt.show()
##
5.497787143782137
m=[[],[],[]] #Liste position pour 3 périodes d'activation
for i in range(0,100): #Boucle ajoute les position pour chaque périodes
m[0].append(euler(0.01,10)[1]) #liste 1 avec periode de 0.01
m[1].append(euler(1,10)[1]) #liste 2 avec periode de 1
m[2].append(euler(10,10)[1])#liste 3 avec periode de 10
##PLOT
fig,ax=plt.subplots(3,2,figsize=(10,7))
ax[0,0].plot(t,euler(0.01,10)[0],'k') # graphique correspondant à liste 1 en fonction du temps
ax[1,0].plot(t,euler(1,10)[0],'k') # graphique correspondant à liste 2 en fonction du temps
ax[2,0].plot(t,euler(10,10)[0],'k') # graphiquecorrespondant à liste 3 en fonction du temps
for i in range(0,3):
ax[i,1].hist(m[i],color='k') # histogramme liste 1 2 3
def plot(i):
ax[i,0].set_xlabel('Temps en périodes') # titre abscisse graphique
ax[i,0].set_ylabel('Position sur x')# titre ordonnée graphique
ax[0,1].set_title('Distribution de la position finale')# titre histogramme
ax[0,0].set_title('Position de la kinésine') # titre graphique
ax[i,1].set_ylabel('Nombre de kinésines')# titre ordonnée histogramme
ax[i,1].set_xlabel('Position finale sur x') # titre abscisse histogramme
ax[i,0].grid()
ax[i,1].grid()
for i in range(0,3): # boucle mise en forme
plot(i)
plt.subplots_adjust(hspace=0.3)
##
l=[[],[],[]] # création de 3 liste pour chaque ligne
for i in range(0,200): # création de 200 points
l[0].append(euler(1,1)[1]) # liste 1 avec amplitude de 1
l[1].append(euler(1,20)[1]) # liste 2 avec amplitude de 10
l[2].append(euler(1,100)[1]) # liste 3 avec amplitude de 100
fig,ax=plt.subplots(3,2,figsize=(10,7))
ax[0,0].plot(t,euler(1,1)[0],'k') #graphique correspondant à liste 1 en fonction du temps
ax[1,0].plot(t,euler(1,20)[0],'k') # graphique correspondant à liste 2 en fonction du temps
ax[2,0].plot(t,euler(1,100)[0],'k') # graphique correspondant à liste 3 en fonction du temps
for i in range(0,3):
ax[i,1].hist(l[i],color='k') # histogramme des listes 1 2 3
def plot(i):
ax[i,0].set_xlabel('Temps en périodes') # titre abscisse graphique
ax[i,0].set_ylabel('Position sur x') # titre ordonnée graphique
ax[0,1].set_title('Distribution de la position finale') # titre histogramme
ax[0,0].set_title('Position de la kinésine') # titre graphique
ax[i,1].set_ylabel('Nombre de kinésines') # titre ordonnée histogramme
ax[i,1].set_xlabel('Position finale sur x') # titre abscisse histogramme
ax[i,0].grid()
ax[i,1].grid()
for i in range(0,3): # boucle mise en forme
plot(i)
plt.subplots_adjust(hspace=0.3)
def eulersym(w,amp): #Resoud l'équation différentielle mais avec un potentiel symétrique
x=0
l=[]
for i in range(0,len(t)):
l.append(x)
x+=(np.random.normal(0,4.5)-vprime(x,amp,1/2)*periode(w)[i])*step
return l,l[-1]
l=[] #Liste vide positions finales
for i in range(0,100): #Boucle ajoute les positions finales à l
l.append(eulersym(1,10)[1])
##PLOT
fig,ax=plt.subplots(1,2,figsize=(15,8))
for i in range(0,5): #Boucle plot les chemins de 10 kinésines
ax[0].plot(t,eulersym(1,10)[0])
ax[0].set_xlabel('Temps en périodes')
ax[0].set_ylabel('Position sur x')
ax[0].set_title('Déplacement de kinésines sur un potentiel symétrique')
ax[1].set_title('Distribution de la position finale')
ax[1].hist(l,color='k') #Histogramme de la distribution de la position finale
ax[1].set_xlabel('Position sur x')
ax[1].set_ylabel('Nombre de kinésines')
ax[1].grid()
ax[0].grid()
plt.show()
##
u = np.linspace(0,20,21) # création d'une liste d'amplitude de 0 à 20
x = [] # creation liste position du premier graphique
for i in range (len(u)): # boucle pour les postion
x.append(euler(1,u[i])[1]/100) # création des position
u1 = np.linspace(0,50,51) # création d'une liste d'amplitude de 0 à 50
x1 = [] # creation liste position du deuxième graphique
for i in range (len(u1)): # boucle pour les postion
x1.append(euler(1,u1[i])[1])# creation liste position final divisé par temps maximal
fig, axs = plt.subplots(2)
axs[0].plot(u, x,'k') # graph 20 amplitude
axs[1].plot(u1, x1,'k') # graphe 50 amplitude
axs[0].set_ylabel('Position sur x') # ordonnées
axs[1].set_xlabel('Temps en périodes')# abscisses
axs[1].set_ylabel('Position sur x') # ordonnées
axs[0].set_title("Position finale en fonction de de l'amplitude") # titre
plt.show
"""
def sorted():
"""
#:return l2 La liste de liste de toutes les position des kinésines à chaque temps
"""
l1 = []
for i in range(0,1000):
l1.append(euler()[0])
l2 = []
for i in range(0, len(l1[0])):
lprov = []
for elem in l1:
lprov.append(elem[i])
l2.append(lprov)
return l2
sort=sorted()
f=open('sort1000.dat','wb') #Crée le fichier sort1000.dat
pickle.dump(sort,f,pickle.HIGHEST_PROTOCOL) #Stocke sort dans sort1000.dat
f.close()
"""
"""
f=open('sort1000.dat','rb') #Ouvre sort1000.dat
d=pickle.load(f) #Lis sort1000.dat
def graph(): #Fonction regroupant les fonction rémanentes dans un plot
ax.set_xlim(-5, 20)
ax.set_ylim(0,1)
plt.grid()
plt.legend()
ax.set_xlabel('Position x')
ax.set_ylabel('Probabilité de présence')
def update_hist(num): #Fonction met à jour l'histogramme
plt.cla()
weight=np.ones_like(d[num])/len(d[0])
plt.axvline(np.mean(d[num]),color='k',linestyle='dashed',label='mean=%1.3e'%np.mean(d[num]))
plt.hist(d[num],ec='k',color='k',weights=weight,label='$\sigma=$%1.3e'%np.std(d[num])+'\nTemps: %1.2f'%(0.01*num),bins=np.arange(min(d[num]), max(d[num]) + 1, 1))
graph()
##PLOT
fig,ax=plt.subplots()
plt.hist(d[0],ec='k',color='k',weights=np.ones_like(d[0])/len(d[0]),label='$\sigma=$%1.3e'%np.std(d[0])+'\nTemps: %1.2f'%(0.01*0),bins=np.arange(min(d[0]), max(d[0]) + 1, 1)) #Affiche l'histogramme au temps 0
graph()
animation = anim.FuncAnimation(fig, update_hist, len(d) ,repeat=False,interval=1) #Fonctio animation
##
matplotlib.rcParams['animation.ffmpeg_path']="C:\\PATHS_\\ffmpeg-2023-04-12-git-1179bb703e-full_build\\bin\\ffmpeg.exe" #Emplacement souhaité du fichier
writervideo=anim.FFMpegWriter(fps=5) #Selectionne le nombre d'images par seconde
animation.save("fokkerplanckfinal1.mp4") #Sauvegarde l'animation
"""
from IPython.display import HTML
# Youtube
HTML('<iframe width="560" height="315" src="https://www.youtube.com/embed/PHOKN2jKSbA?rel=0&controls=0&showinfo=0" frameborder="0"allowfullscreen></iframe>')
/shared-libs/python3.9/py-core/lib/python3.9/site-packages/IPython/core/display.py:419: UserWarning: Consider using IPython.display.IFrame instead
warnings.warn("Consider using IPython.display.IFrame instead")
f=open('sort1000.dat','rb')
d=pickle.load(f)
print(d)
std=[]
mean=[]
for i in range(0,len(d)):
std.append(np.std(d[i]))
mean.append(np.mean(d[i]))
t=np.arange(0,10000,1)
def modsqrt(t,D):
return np.sqrt(2*D*t)
def modlin(t,a,b):
return a*t+b
D=curve_fit(modsqrt,t,std)[0]
a,b=curve_fit(modlin,t,mean)[0]
def recurplot(i):
ax[i].legend()
ax[i].grid()
fig,ax=plt.subplots(1,2,figsize=(10,5))
ax[0].plot(mean,'ok',label='Données',markersize=1)
ax[0].plot(t,modlin(t,a,b),label='Modèle linéaire: a=%2.3e'%a)
ax[1].plot(std,'ok',label='Données',markersize=1)
ax[1].plot(t,modsqrt(t,D),'r',label='Modèle sqrt: D=%2.3e'%D)
ax[0].set_title('Position moyenne au cours du temps')
ax[1].set_title('Ecart-type au cours du temps')
ax[0].set_xlabel('Temps')
ax[0].set_ylabel('Position sur x')
ax[1].set_ylabel('Position sur x')
ax[1].set_xlabel('Temps')
for i in range(0,2):
recurplot(i)
plt.show()
IOPub data rate exceeded.
The notebook server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--NotebookApp.iopub_data_rate_limit`.
Current values:
NotebookApp.iopub_data_rate_limit=1000000.0 (bytes/sec)
NotebookApp.rate_limit_window=3.0 (secs)
/tmp/ipykernel_93/3608389216.py:13: RuntimeWarning: invalid value encountered in sqrt
return np.sqrt(2*D*t)
https://anatomypubs.onlinelibrary.wiley.com/doi/full/10.1002/ar.20599 Cellular Motors for Molecular Manufacturing (2007) C.Z. Dinu, D.B. Chrisey, S. Diez, J. Howard
https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=7745727529a1c5c4e2bf91f2171361e2dce13cc3 FORCE GENERATION BY CELLULAR MOTORS (2003) FRIEDRICH WANKA and EVERARDUS J.J. VAN ZOELEN
https://link-springer-com.accesdistant.sorbonne-universite.fr/article/10.1385/CBB%3A38%3A2%3A191 Brownian ratchet models of molecular motors(2003) Rachid Ait-Haddou & Walter Herzog
https://link-springer-com.accesdistant.sorbonne-universite.fr/content/pdf/10.1007/978-3-030-72515-0_4?pdf=chapter%20toc Stochastic Processes in Cell Biology (2022) Paul C. Bresslof