The goal of this project was to modellise the polarisation of /C. elegans/ once-cell embryo. In the following we construct step by step a simple model and explore some of its extensions.
Overview:
We will use a discretisation of time and space. The concentration of each species will be described by an EDP.
In [1]:
# Some boilerplate code
from __future__ import division
import numpy as np
import pandas
import matplotlib.pyplot as plt
%matplotlib inline
%cd data/
cytosol,advection = None,None
plt.rcParams['figure.figsize'] = 15,7
These are the measured physical constants
In [2]:
def get_P_param(**kargs):
param = {}
#Discretisation
param["L"] = 134.6 #µm
param["T"] = 700 #s
# Kinetics
param["kon"] = 4.74e-2 # μm s−1
param["koff"] = 7.3e-3 # s-1
param["psi"] = 0.174 #µm-1
# Concentration
param["rho"] = 1 #µm--3
# Diffusion
param["D"] = 0.15 #µm^2.s-1
# Advection parameter
param["xi"] = 1 #adimensional
#interaction
param["alpha"] = 1
param["beta"] = 2
param["kap"] = 0.19
param["kpa"] = 0.3
for k,v in kargs.items():
param[k] = v
return param
def get_A_param(**kargs):
param = get_P_param()
param["kon"] = 8.58e-3 # µm s-1
param["koff"] = 5.4e-3 # s-1
param["xi"] = 1 # adimensional
param["rho"] = 1.56 #µm--3
param["D"] = 0.28
for k,v in kargs.items():
param[k] = v
return param
# If you want to get a new set of parameters just call:
theta = get_P_param()
# You can specifie the value of some of the parameters:
theta = get_A_param(a=4,D=1000)
Model implementation
In [3]:
def play(p,param,components):
""" The simplest form of the model: Diffusion only
Args:
p (np.array): A timexSpace array, the first line give
the initial condtion, its shape gives the number of
discretization volumes.
param (dict): Parameters
components (tuple): the list of functions to use in the computation of dP.
Returns:
(np.array): with a timestep by line.
"""
dt = param["T"]/float(p.shape[0]-1) #Time discretisation.
l = param["L"]/float(p.shape[1]-1) #Space discretisation.
param["dt"] = dt
param["l"] = l
print "Simulation of {}, T={} (dt={}), L={} (dx={})".format([x.__name__ for x in components],param["T"],dt,param["L"],l)
#Cytosol
param["cyto"] = np.zeros(p.shape[0])
param["cyto"][0] = param["rho"] - param["psi"]/p.shape[1] * np.sum(p[0,:])
# Checking for advection
if advection in components:
dmax = param["F"].max() * param["dt"]
if dmax >= 0.1*l:
print ("Warning ! Maximal flow movment in one step is {}µm "
"whereas the space discretization length is {}µm.\n "
"You could try to reduce dt to improve that.").format(dmax,l)
# Time loop
for t in range(p.shape[0]-1):
dpoverdt = []
for fct in components:
dpoverdt.append(fct(p[t,:],param,dt,l,t))
p[t+1,:] = (p[t,:] + dt*(np.sum(dpoverdt,0)))
if cytosol in components:
param["cyto"][t+1] = param["cyto"][t] - dt*np.sum(cytosol(p[t,:],param,dt,l,t))
return p,param
In [4]:
def display(p,param={},cyto=False):
"""Display function
Args:
p (np.array): A timeXspace array.
param (dict): the parameters
cyto (bool): display the membrane/cytosol concentration"""
plt.subplot(1,2,1)
plt.imshow(p,interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0))
plt.xlabel("distance (um)")
plt.ylabel("time (sec)")
plt.subplot(1,2,2)
plt.xlabel(r'Distance ($\mu m$)')
plt.ylabel("Concentration ($.m^{-3}$)")
plt.plot([x*param["l"]-param["L"]/2 for x in range(p.shape[1])],p[-1,:],label="t=final")
plt.plot([x*param["l"]-param["L"]/2 for x in range(p.shape[1])],p[0,:],label="t=0")
plt.plot([x*param["l"]-param["L"]/2 for x in range(p.shape[1])],p[1,:],label="t=1")
plt.legend()
if cyto:
plt.show()
j=200
plt.plot(param["cyto"][:j],label="Cytosol")
plt.plot(p.sum(1)[:j],label="Membrane")
plt.plot((p.sum(1)+param["cyto"])[:j],label="Overall")
plt.xlabel("time (sec)")
plt.ylabel("$P$")
plt.legend()
In the following you have the simplest implementaton of the model, using only diffusion.
$$\frac{\partial P_i}{\partial t} = -D\nabla P_i $$Discretised by looking at each eulerian control volume:
$$ P_i(t+dt) = P_i(t) + dt \frac{D}{l^2} \left(P_{i-1}(t) - 2P_i(t) + P_{i+1}\right)$$
In [5]:
def diffusion(p,param,dt,l,t):
return (param["D"] * l**-2 *
(np.concatenate((p[1:],[p[0]]))
- 2*p[:]
+ (np.concatenate(([p[-1]],p[:-1])))))
In [191]:
# First we definie an empy matrix (time*space) that we will fill step by step.
p = np.zeros((3000, 65))
# Initial conditions : a dirac in the middle.
p[0,int(p.shape[1]/2)] = 1
# And two others points because why not.
p[0,5] = 0.5
p[0,50:60] = 0.3
# Run the simulation !
p,param = play(p,get_A_param(T=300),[diffusion])
display(p,param)
plt.savefig("03_diffusion.svg")
In [7]:
# The same thing but starting from a sinusoidal wave.
p = np.zeros((3000, 65))
p[0,:] = 1+np.sin(np.arange(p.shape[1]))
# Run the simulation !
p,param = play(p,get_P_param(T=100),[diffusion])
display(p,param)
Convergence plots are of critical importance in order to assess the correctness of our discretization scheme and choose the time and space step size.
The idea is to consider a simple case where an analytical solution is known (for instance the diffusion of a dirac distribution gives a gaussian) and compare it to the result given by our model for different resolution.
$$G(x,t) = \frac{1}{\sqrt{4 \pi D t}} exp \left(-\frac{(x-x_0)^2}{4Dt}\right)$$Our space discretisation scheme should lead to an error of $O(dx^2)$. We can prove it using a Taylor expansion of the concentration around the focal point.
$$P(x+dx) = P(x) + \frac{dx}{1!} \frac{\partial P}{\partial x}\bigg |_x + \frac{dx^2}{2!} \frac{\partial^2 P}{\partial x^2}\bigg |_x + \frac{d x^3}{3!} \frac{\partial^3 P}{\partial x^3}\bigg |_x$$$$P(x-dx) = P(x) - \frac{d x}{1!} \frac{\partial P}{\partial x}\bigg |_x + \frac{dx^2}{2!} \frac{\partial^2 P}{\partial x^2}\bigg |_x - \frac{d x^3}{3!} \frac{\partial^3 P}{\partial x^3}\bigg |_x$$
By a linear combination (1)+(2) we get:
$$ \frac{\partial^2 P}{\partial x^2}\bigg |_x = \frac{P(x+dx) + P(x-dx)-2P(x)}{d x^2} + \frac{1}{d x^2}\frac{2d x^4}{4!} \frac{\partial^4 P}{\partial x^4}$$$$ \frac{\partial^2 P}{\partial x^2}\bigg |_x = \frac{P(x+dx) + P(x-dx)-2P(x)}{d x^2} +O(dx^2)$$Conversely, our time discretisation scheme should lead to an error of $O(dt)$.
(1) leads to: $$P(t+dt) = P(t) + dt\frac{\partial P}{\partial t}\bigg |_t + O(dt^2)$$ $$\frac{P(t+dt) - P(t)}{dt} = \frac{\partial P}{\partial t}\bigg |_t + O(dt)$$
In [8]:
def analytical_vector(N,param,disp=True):
"""
Analatical solution of the dirac's diffusion
Args:
N (int): number of grid points
param (dict): A dictionnary of parameters with at least
"D", the diffusion coefficient and "T" the time.
"""
a = (4*np.pi*param["D"]*param["T"])**(-0.5)
b = 4*param["D"]*param["T"]
x_0 = int(N/2)
c = np.array([-(x-x_0)**2 for x in np.arange(N)])
out = a * np.exp(c/b)
if disp:
print "Analytical results, T={}, L={} (dx={})".format(param["T"],param["L"],param["L"]/N)
return out
In [187]:
#######################
#Convergence in SPACE...
#######################
e = [] #error
e2 = [] #error
Nvector = [] #number of steps
runtime = 10
plt.subplot(1,2,1)
for i,N in enumerate([int(10**x) for x in np.arange(1,3,0.2)]):
#for i,N in enumerate([500]):
# simulation
p = np.zeros((3000,N))
p[0,:] = analytical_vector(N,get_P_param(T=10),False)
p,param = play(p,get_P_param(T=runtime),[diffusion])
# analytical solution
analytical = analytical_vector(N,get_P_param(T=10+runtime))
# save the results
Nvector.append(p.shape[1])
e.append(np.max(abs(p[-1,:]-analytical)))
# plots
if i%2==0:
plt.plot( [x*param["l"]-param["L"]/2 for x in range(p.shape[1])], analytical-p[-1,:], alpha=0.5,
label="{}".format(N))
plt.legend(title="Gridpoints")
plt.xlabel(r'Distance ($\mu m$)')
plt.ylabel("Error")
# display
plt.subplot(1,2,2)
plt.plot(Nvector,e,label="Tfinal")
#plt.plot(Nvector,e2,label="Tbegining")
plt.plot(Nvector[:6],[100*x**-2 for x in Nvector[:6]],label="$y=100x^{-2}$")
CFL = param["L"]/np.sqrt(2*param["D"]*param["dt"])
plt.legend()
plt.xlabel("Number of grid points")
plt.ylabel("Maximal Error")
plt.yscale('log')
plt.xscale('log')
plt.vlines([64,150,CFL],*plt.ylim(),label="CFL limit")
ax = plt.gca()
ax.annotate(r'$CFL: \frac{L}{\sqrt{2Ddt}}$', xy=(CFL,0.2),xycoords="data", textcoords="offset points", xytext=(-80,0),
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
ax.annotate(r'Data', xy=(64,0.3),xycoords="data", textcoords="offset points", xytext=(-40,0),
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
ax.annotate('Time convergence plot', xy=(150,0.3),xycoords="data", textcoords="offset points", xytext=(+50,+10),
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.savefig("04_convergence_space.svg")
plt.show()
In [14]:
def diff_in_time(N):
a = [analytical_vector(N,get_P_param(T=x),False)[N/2] for x in np.arange(1,100,1)]
p = np.zeros((100,N))
p[0,N/2]=1
p,param = play(p,get_P_param(T=100),[diffusion])
x = p[1:,N/2]
plt.plot(range(len(a)),a,label="analytical")
plt.plot(range(len(x)),x,label="experimental")
plt.xlabel("Time")
plt.ylabel("P(x=0)")
#plt.yscale("log")
plt.ylim((-1,1))
plt.legend(title="N={}".format(N),loc="lower left")
plt.subplot(1,3,1)
diff_in_time(10)
plt.subplot(1,3,2)
diff_in_time(150)
plt.subplot(1,3,3)
diff_in_time(250)
In [186]:
#######################
#Convergence in TIME...
#######################
#Analytical solution
e = [] #error
Nvector = [] #number of steps
Nx = 150
plt.subplot(1,2,1)
for i,N in enumerate([int(10**x) for x in np.arange(1,3.8,0.2)]):
# analytical solution
analytical = analytical_vector(Nx,get_P_param(T=20))
# simulation
p = np.zeros((N,Nx))
p[0,:] = analytical_vector(Nx,get_P_param(T=10),False)
p,param = play(p,get_P_param(T=10),[diffusion])
# Compute the error and save the results for plotting.
Nvector.append(p.shape[0])
e.append(np.max(abs(p[-1,:]-analytical)))
# plots
plt.plot( [x*param["l"]-param["L"]/2 for x in range(p.shape[1])], analytical-p[-1,:], alpha=0.5,
label="{}".format(N))
plt.legend(title="Timesteps")
plt.xlabel("Space")
plt.xlabel(r'Distance ($\mu m$)')
# display
plt.subplot(1,2,2)
plt.plot(Nvector,e)
#plt.legend()
plt.xlabel("Number of timesteps")
plt.ylabel("Maximal Error")
plt.vlines([3000],*plt.ylim())
ax = plt.gca()
ax.annotate(r'Space convergence plot', xy=(3000,0.008),
xycoords="data", textcoords="offset points", xytext=(-200,0),
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.savefig("04_convergence_time.svg")
plt.show()
We add the cytosol compartment and the reaction of attachment/falling of the membrane.
$$P \xrightarrow{k_{off}} P_{cyto}$$$$P_{cyto} \xrightarrow{k_{on}} P$$The concentration of particles in the cytosol can be computed from the total concentration of proteins in the cell: $$\rho_P = P_{cyto} + P_{membrane}$$ $$\rho_P = P_{cyto} + \frac{\psi}{L} \int_L l \cdot P dx$$ $$\rho_P = P_{cyto} + \frac{\psi \cdot l}{L} \int_L P dx$$ $$\rho_P = P_{cyto} + \frac{\psi}{N} \int_L P dx$$
If we discretize this over our control volumes :
$$\rho_P = P_{cyto} + \frac{\psi}{N} \sum P_i \Leftrightarrow P_{cyto} = \rho_P - \frac{\psi}{N} \sum P_i $$So the model is now:
$$ P_i(t+dt) = P_i(t) + dt [\frac{D}{l^2} \left(P_{i-1}(t) - 2P_i(t) + P_{i+1}\right) + k_{on}P_{cyto}\frac{\psi}{N} - k_{off}P_i(t)] $$$$ P_{cyto}(t+dt) = P_{cyto}(t) + dt\left(-k_{on}P_{cyto}(t)\frac{\psi}{N} +k_{off}P_{membrane}(t)\right) $$
In [17]:
def cytosol(p,param,dt,l,t):
return param["cyto"][t]*param["kon"]*param["psi"] - p*param["koff"]
In [7]:
# First we definie an empy matrix (time*space) that we will fill step by step.
p = np.zeros((3000, 65))
# Run the simulation !
p,param = play(p,get_P_param(rho=1),[cytosol])
display(p,param,True)
plt.savefig("03_turnover.svg")
In [44]:
p = np.zeros((3000, 65))
p[0,int(p.shape[1]/2)] = 1 #A dirac
# Run the simulation !
#p,param = play(p,get_P_param(rho=0),[cytosol])
p,param = play(p,get_P_param(rho=0),[diffusion,cytosol])
display(p,param,True)
The expression of the temporal variation of concentration due to the advection is of the form:
$$ \frac{\partial P}{\partial t} = \nabla \cdot \vec{f}P $$As the flow is compressible, $\nabla \cdot \vec f \neq 0$, so we cannot simplify further this expression.
When we discretize our model is now:
$$ P_i(t+dt) = P_i(t) + dt [\frac{D}{l^2} \left(P_{i-1}(t) - 2P_i(t) + P_{i+1}\right) + k_{on}P_{cyto}\frac{\psi}{N} - k_{off}P_i(t)] + -F[] $$$$ P_{cyto}(t+dt) = P_{cyto}(t) + dt\left(-k_{on}P_{cyto}(t)\frac{\psi}{N} +k_{off}P_{membrane}(t)\right) $$
In [16]:
def advection(p,param,dt,l,t):
return ( np.concatenate(([param["F"][t,-1]],param["F"][t,:-1])) * l**-1 * param["xi"] * np.concatenate(([p[-1]],p[:-1]))
- np.concatenate((param["F"][t,1:],[param["F"][t,0]])) * l**-1 * param["xi"] * np.concatenate((p[1:],[p[0]]))
)
In [1]:
def artificial_flow(Nt,Nx,s=1e-2):
flow_artificial = np.zeros((Nt, Nx)) - s
flow_artificial[:,int(Nx/2):] = -flow_artificial[:,int(Nx/2):]
flow_artificial[:,:int(0.2*Nx)] = 0
flow_artificial[:,-int(0.2*Nx):] = 0
flow_artificial[:,int(Nx/2)-0.05*int(Nx/2):int(Nx/2)+0.05*int(Nx/2)] = 0
return flow_artificial
#plt.subplot(1,2,1)
#plt.imshow(artificial_flow(2,150),interpolation="none",
# aspect="auto",cmap="rainbow",
# extent=(-0.5*param["L"],0.5*param["L"],param["T"],0))
#plt.colorbar()
#plt.xlabel("distance")
#plt.xlabel("time")
#plt.savefig("09_flow_cst.png")
In [215]:
plt.subplot(1,2,1)
plt.imshow(artificial_flow(2,150),interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0))
plt.subplot(1,2,2)
plt.imshow(flow,interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0))
plt.plot()
Out[215]:
In [2]:
p = np.zeros((9000, 150))
p[0,:] = analytical_vector(150,get_P_param(T=20),False)
p,param = play(p,get_P_param(F=artificial_flow(*p.shape,s=1e-4)),[advection])
display(p,param)
In [194]:
# Get the experimental values (these are in µm.s-1 so we multiply by 60 to have seconds)
p = np.zeros((flow.shape[0], 65)) + 1/65
#p[0,:] = analytical_vector(65,get_P_param(T=10),False)
p,param = play(p,get_P_param(F=flow,rho=0),[advection])
display(p,param)
plt.savefig("03_advection.svg")
In [11]:
def play2(p,a,param,param_a,components,components2=[]):
""" The simplest form of the model: Diffusion only
Args:
p (np.array): A timexSpace array, the first line give
the initial condtion, its shape gives the number of
discretization volumes.
a (np.array): A timexSpace array, the first line give
the initial condtion, its shape gives the number of
discretization volumes.
param (dict): Parameters for p
param_a (dict): Parameters for a
components (tuple): the list of functions to use in the computation of dP/dA.
components_2 (tuple): the list of functions to use in the computation of dP&dA.
Returns:
(np.array): with a timestep by line.
"""
dt = param["T"]/float(p.shape[0]-1) #Time discretisation.
l = param["L"]/float(p.shape[1]-1) #Space discretisation.
param["dt"] = dt
param["l"] = l
print "Simulation of 2 species with {},{}, T={} (dt={}), L={} (dx={})".format([x.__name__ for x in components],
[x.__name__ for x in components2],
param["T"],dt,param["L"],l)
#Cytosol
param["cyto"],param_a["cyto"] = np.zeros(p.shape[0]), np.zeros(a.shape[0])
param["cyto"][0] = param["rho"] - param["psi"]/p.shape[1] * np.sum(p[0,:])
param_a["cyto"][0] = param_a["rho"] - param_a["psi"]/a.shape[1] * np.sum(a[0,:])
# Checking for advection stability
if advection in components:
dmax = max(param["F"].max() * param["dt"],param_a["F"].max() * param["dt"])
if dmax >= 0.1*l:
print ("Warning ! Maximal flow movment in one step is {}µm "
"whereas the space discretization length is {}µm.\n "
"You could try to reduce dt to improve that.").format(dmax,l)
# Time loop
for t in range(p.shape[0]-1):
dpoverdt = []
daoverdt = []
for fct in components:
dpoverdt.append(fct(p[t,:],param,dt,l,t))
daoverdt.append(fct(a[t,:],param_a,dt,l,t))
if cytosol in components:
param["cyto"][t+1] = param["cyto"][t] - dt*np.sum(cytosol(p[t,:],param,dt,l,t))
param_a["cyto"][t+1] = param_a["cyto"][t] - dt*np.sum(cytosol(a[t,:],param_a,dt,l,t))
for fct in components2:
x,y,param,param_a = fct(p[t,:],a[t,:],param,param_a,dt,l,t)
dpoverdt.append(x)
daoverdt.append(y)
p[t+1,:] = (p[t,:] + dt*(np.sum(dpoverdt,0)))
a[t+1,:] = (a[t,:] + dt*(np.sum(daoverdt,0)))
return p,a,param,param_a
In [10]:
def display2(p,a,param,param_a,cyto=False,flow=False):
"""Display function
Args:
p (np.array): A timeXspace array.
param (dict): the parameters
cyto (bool): display the membrane/cytosol concentration"""
vmin = min(p.min(),a.min())
vmax = max(p.max(),a.max())
plt.subplot(1,2,1)
plt.imshow(p,interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0),
vmin=vmin,vmax=vmax)
plt.title("P")
plt.xlabel("distance (um)")
plt.ylabel("time (sec)")
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(a,interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0),
vmin=vmin,vmax=vmax)
plt.xlabel("distance (um)")
plt.ylabel("time (sec)")
plt.title("A")
plt.colorbar()
if cyto:
plt.show()
plt.subplot(2,1,1)
plt.plot(param["cyto"],label="P Cytosol")
plt.plot(p.sum(1),label="P Membrane")
plt.plot(p.sum(1)+param["cyto"],label="P Overall")
print max(p.sum(1)+param["cyto"]),min(p.sum(1)+param["cyto"])
plt.xlabel("time (sec)")
plt.ylabel("Concentration")
plt.legend()
plt.subplot(2,1,2)
plt.plot(param_a["cyto"],label="A Cytosol")
plt.plot(a.sum(1),label="A Membrane")
plt.plot(a.sum(1)+param_a["cyto"],label="A Overall")
plt.xlabel("time (sec)")
plt.ylabel("Concentration")
plt.legend()
if flow:
plt.show()
plt.subplot(1,2,2)
plt.imshow(param["F"],interpolation="none",
aspect="auto",cmap="rainbow",
extent=(-0.5*param["L"],0.5*param["L"],param["T"],0),
vmin=vmin,vmax=vmax)
plt.xlabel("distance (um)")
plt.ylabel("time (sec)")
plt.title("Flowspeed um.sec-1")
For which values of $\xi$ would an embryo survive ?
4 conditions for survival:
We define a survival condition as $\forall t \in [500-700]$:
So we have a survivability index as the product of the squared errors of those vecors normalised by the one of the data.
$$E_{C,\theta} = \left( \frac{(SE_{model}(C,\theta)-SE_{data}(C))^2 }{SE_{data}(C)} \right)^2$$with $SE(C,\theta)$ the standard error for the condition $C$ and model parameters $\theta$, The survivability index is $S(\theta)$
$$S(\theta) = \frac{1}{200} \sum_{t = 500}^{700-dt} E_{Stability,\theta} \cdot E_{Polarity,\theta} \cdot E_{Anticorrelation,\theta} $$If the elements we defined are indeen sufficient conditions for the embryo to survive, it should survive if $S(\xi) \approx 1$ and die if $S(\xi) \approx 0$
In [31]:
def error(model,param,data,t0=500):
S = 0
Epola = 0
Eantico = 0
Estab = 0
dtd=1
t_range_data = np.arange(50,70,dtd)
i = param["dt"]/500
l = model["a"].shape[0]-i
dtm = int(l/20)
t_range_model = np.arange(i,model["a"].shape[0],dtm)
for tm,td in zip(t_range_model,t_range_data):
#antico
se_data = np.sum(np.square(data["a"][td,:]-data["p"][td,:]))
se_model = np.sum(np.square(model["a"][tm,:]-model["p"][tm,:]))
e_antico = ((se_model - se_data)/(se_data))**2
Eantico += e_antico
#stability
se_data = np.sum(np.square(data["p"][td+dtd,:]-data["p"][td,:]))
se_model = np.sum(np.square(model["p"][tm+dtm,:]-model["p"][tm,:]))
e_stab = ((se_model - se_data)/(se_data))**2
Estab += e_stab
#polarity
se_data = np.sum(np.square((np.concatenate((data["p"][td,:int(data["p"].shape[1]/2)+1:-1],data["p"][td,-int(data["p"].shape[1]/2)::-1]))
-data["p"][td,:])))
se_model = np.sum(np.square((np.concatenate((model["p"][td,:int(model["p"].shape[1]/2)+1:-1],model["p"][td,-int(model["p"].shape[1]/2)::-1]))
-model["p"][td,:])))
e_pola = ((se_model - se_data)/(se_data))**2
Epola += e_pola
S += e_antico * e_pola * e_stab
return S/len(t_range_model),Epola/len(t_range_model),Eantico/len(t_range_model),Estab/len(t_range_model)
In [33]:
e = []
Epola = []
Eantico = []
Estab = []
xi_range = np.arange(0,1,0.25)
for xi in xi_range:
p = np.zeros(flow.shape)
a = np.zeros(flow.shape)
p,a,param,param_a = play2(p,a,get_P_param(F=flow,xi=xi),get_A_param(F=flow),[diffusion,cytosol,advection])
x,y,z,t = error({"a":a,"p":p},param,{"a":apar,"p":ppar},t0=500)
e.append(x)
Epola.append(y)
Eantico.append(z)
Estab.append(t)
#plt.plot(xi_range,xi_range)
plt.plot(xi_range,e)
plt.plot(xi_range,Epola,label="Polarity")
plt.plot(xi_range,Eantico,label="Anticorrelation")
plt.plot(xi_range,Estab,label="Stability")
plt.xlabel("xi_P")
plt.ylabel("S(xi_P)")
plt.legend()
Out[33]:
In [59]:
def interaction(p,a,param,param_a,dt,l,t):
dp = - param["kpa"] * p * a**param["alpha"]
da = - param_a["kap"] * p**param_a["beta"] * a
param["cyto"][t+1] = param["cyto"][t] - dt*np.sum(dp)
param_a["cyto"][t+1] = param_a["cyto"][t] - dt*np.sum(da)
return dp,da,param,param_a
In [182]:
p = np.zeros((3000,150))+1/150
a = np.zeros((3000,150))+1.56/150
display2(*play2(p,a,get_P_param(kpa=0.2,xi=1,rho=0),get_A_param(F=flow,rho=0),[cytosol],[interaction]),cyto=True)
In [229]:
# First we definie an empy matrix (time*space) that we will fill step by step.
p = np.zeros(flow.shape)
a = np.zeros(flow.shape)
# Run the simulation !
## Curated public version of the notebook: upublished data are not yet avaiable.
#flow = pandas.read_csv("flow-myosin-um-per-min_interpolated_40.txt").values / 60
display2(*play2(p,a,get_P_param(F=artificial_flow(*flow.shape,s=5e-2),kpa=0.1,xi=1),get_A_param(F=artificial_flow(*flow.shape,s=5e-2)),[diffusion,cytosol,advection],[]))
plt.savefig("09_cst_flow.png")
plt.show()
p = np.zeros(flow.shape)
a = np.zeros(flow.shape)
# ## Curated public version of the notebook: upublished data are not yet avaiable.
# flow = pandas.read_csv("flow-myosin-um-per-min_interpolated_40.txt").values / 60
display2(*play2(p,a,get_P_param(F=flow,kpa=0.1,xi=1),get_A_param(F=flow),[diffusion,cytosol,advection],[]))
plt.savefig("05_onlyflow.png")
plt.show()
## Curated public version of the notebook: upublished data are not yet avaiable.
#print "Experimental results:"
#display2(ppar,apar,get_P_param(),get_A_param())
#plt.savefig("05_data.png")