%load_ext autoreload
%autoreload 2
%matplotlib inline
%run nb_init.py
2015-01-26 09:33:26,838 -leg_joint -INFO -successfully imported leg_joint
According to the Farhadifar et al. paper, in the case of a regular hexagonal latice, energy is given by:
$$ E = N\frac{K}{2} (A - A_0)^2 + N\frac{\Gamma}{2} L^2 + 6N\frac{\Lambda}{2}\ell $$In our model, the area dependant term is replaced by a volume term : $$ E = N\frac{K_v}{2} (V - V_0)^2 + N\frac{\Gamma}{2} L^2 + 3N\Lambda\ell $$
Here the $V = A(\rho - R_L) = Ah $ where $R_L$ is the lumen radius and $\rho$ the radius of the apical sheet of the epithelium. As they did, we define the adimentional contractility $\bar\Gamma = \Gamma/K_vA_0h_0^2$ and line tension $\bar\Lambda = \Lambda /K_v (A_0^{3/2}h_0^2)$, where $h_0$ is such that $V_0 = A_0h_0$.
As the epithelium is distributed over a cylinder, the radius of the cylinder determines the number of cells forming a ring around it, $ N_\perp $, so we have: $ \ell = \frac{2\pi}{N_\perp}\rho $ if the cells are aligned along the $z$ axis. The perimeter $L$ of a cell is equal to $6\ell$ and the area $A$ equals $(3\sqrt{3}/2)\ell^2$. We define the constant $\mu = 6\left(2/3\sqrt{3}\right)^{1/2}$, then $\ell = \frac{\mu}{6}A^{1/2}$ and $L^2 = \mu^2 A$.
The normalized energy $\bar E = E/NK_v(A_0h_0)^2$ then reads: $$ \bar E = \frac{1}{2} \left(\frac{V - V_0}{V_0}\right)^2
\bar E = \frac{1}{2} \left(\delta^3 - 1\right)^2
\bar E_{2D} = \frac{1}{2} (\delta^2 - 1)^2
$$
We can implement those functions rightaway:
import numpy as np
import matplotlib.pylab as plt
mu = 6 * np.sqrt(2. / (3 * np.sqrt(3)))
def elasticity(delta):
return (delta**3 - 1 )**2 / 2.
def contractility(delta, gamma):
return gamma * mu**2 * delta**2 / 2.
def tension(delta, lbda):
return lbda * mu * delta / 2.
def isotropic_energy(delta, gamma, lbda):
energy = (elasticity(delta)
+ contractility(delta, gamma)
+ tension(delta, lbda))
return energy
def isotropic_energy_2D(delta, gamma, lbda):
elasticity = (delta**2 - 1)**2 / 2.
contractility = gamma * mu**2 * delta**2 / 2.
tension = lbda * mu * delta / 2.
return elasticity + contractility + tension
delta = np.linspace(0, 1.2, 200)
#aspects = np.logspace(-6, 0, 40, base=2)
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(delta, isotropic_energy_2D(delta, 0.04, 0.12), 'b-', lw=2, alpha=0.8, label='2D case')
ax.plot(delta, isotropic_energy(delta, 0.04, 0.12), 'k-', lw=2, alpha=0.8, label='3D case')
#for aspect in aspects:
# ax.plot(delta, isotropic_energy(delta, aspect, 0.04, 0.12), 'b-',
# label=(r'$R_L/R_0 = %.3e$' % aspect), alpha= aspect * 0.6)
ylbl = ax.set_ylabel(r'Normalized energy $\bar E$')
xlblb = ax.set_xlabel('')
xlblb = ax.set_xlabel(r'Scaling factor $\delta$')
ax.legend(loc='upper left')
fig.savefig(lj.data.get_image('isotropic_energy_vs_scaling.svg'))
/home/guillaume/anaconda/envs/python3/lib/python3.4/site-packages/matplotlib-1.4.x-py3.4-linux-x86_64.egg/matplotlib/figure.py:1644: UserWarning: This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect. warnings.warn("This figure includes Axes that are not "
The ground state of the regular epithelium is attained when enregy is minimal with respect to $\delta$:
$$ \frac{\partial \bar E}{\partial \delta} = 0 $$That is (dropping the top bar notation):
$$ 3 \delta^2(\delta^3 - 1) + \mu^2\Gamma\delta + \frac{\mu\Lambda}{2} = 0 $$In 2D:
$$ 2 \delta(\delta^2 - 1) + \mu^2\Gamma\delta + \frac{\mu\Lambda}{2} = 0 $$Once again let's code this:
def isotropic_grad_poly(lbda, gamma):
grad_poly = [3, 0, 0,
-3,
mu**2 * gamma,
mu * lbda / 2.]
return grad_poly
def isotropic_grad(delta, lbda, gamma):
grad_poly = isotropic_grad_poly(lbda, gamma)
return np.polyval(grad_poly, delta)
def isotropic_grad_poly_2D(lbda, gamma):
grad_poly = [2, 0,
mu**2 * gamma - 2,
mu * lbda / 2.]
return grad_poly
def isotropic_grad_2D(delta, lbda, gamma):
grad_poly = isotropic_grad_poly_2D(lbda, gamma)
return np.polyval(grad_poly, delta)
## The (archetypal) case I simulations in the article are for
## the following values:
gamma_case1 = 0.04
lbda_case1 = 0.12
# We can plot the above function for a range of values for aspect
deltas = np.linspace(0, 1.1, 200)
fig, ax = plt.subplots(1, 1, figsize=(6.,4.))
grad = isotropic_grad(deltas, lbda_case1, gamma_case1)
line = ax.plot(deltas, grad, 'k-', alpha=0.8, lw=2, label='3D case')
grad2D = isotropic_grad_2D(deltas, lbda_case1, gamma_case1)
line = ax.plot(deltas, grad2D, 'b-', alpha=0.8, lw=2, label='2D case')
zero_x = ax.plot(deltas, np.zeros_like(deltas), 'k-')
ylbl = ax.set_ylabel(r'Isotropic gradient amplitude'
'$\partial E / \partial \delta$')
xlblb = ax.set_xlabel(r'Scaling factor $\delta$')
ax.legend(loc='upper left')
plt.draw()
fig.savefig(lj.data.get_image('isotropic_gradient_vs_scaling.svg'))
In order for the apical junction to find a stable configuration, the energy must have a minimum as a function of $\delta$. We can look at several cases and deduce the necessary conditions on $\bar\Gamma$ and $\bar\Lambda$ for this to apply.
The tissue will behave as a soft network if the energy minimum for equilibrium volumes equal the initial volume ($\delta = 1$), and contractility and elasticity can compensate each over. There exist a minimum energy in this case if: $$ \mu^2 \Gamma \delta + \mu\Lambda / 2 < 0 $$ The boundary between soft and cristal-like network is given by the line $\Gamma = - 2 \Lambda / 2\mu$ in the $(\Gamma, \Lambda)$ plane.
To fix the correct conditions we have to look in the $(\Gamma, \Lambda)$ plane for real values of the gradient roots.
We can seek the maximum possible values of $\Gamma$ for a given $\Lambda$. We have: $$ \Gamma = - \frac{\Lambda}{2\mu\delta} + \frac{3}{\mu^2}\delta(1 - \delta^3) $$
$\Gamma$ is maximal with respect to $\delta$ if $$ \begin{align*} \frac{\partial \Gamma}{\partial \delta} = \frac{\Lambda}{2\mu\delta^2} + 3\frac{1 - 4\delta^3}{\mu^2} = 0\\ \Leftrightarrow -12\delta^5 + 3\delta^2 + \Lambda\mu / 2 = 0 \end{align*} $$
The phase space boundary for valid values of $\Gamma$ and $\Lambda$ is given by the values of $\Gamma$ for the roots of the above polynomial as a function of $\Lambda$.
If $\Lambda = 0$, this gives $\delta_m(\Lambda = 0) = 4^{-1/3}$ and $\Gamma_m(\Lambda = 0) \simeq 0.123$. With a similar reasoning for $\Lambda$, we have $\delta_m(\Gamma = 0) = (2 / 5)^{1/3}$ and $\Lambda_m(\Gamma = 0) \simeq 0.525$
We can solve numerically for the other values of $\Lambda$
### This function does a brute force search for the roots of
## the gradient in the lambda, gamma plane
def find_grad_roots(lbda, gamma):
p = isotropic_grad_poly(lbda, gamma)
roots = np.roots(p)
good_roots = np.real([r for r in roots if np.abs(r) == r])
np.sort(good_roots)
if len(good_roots) == 1:
return good_roots
elif len(good_roots) > 1:
return good_roots[0]
else:
return np.nan
def find_boundary_roots(lbda):
delta_poly = [-12, 0, 0, 3, 0, lbda * mu / 2]
roots = np.roots(delta_poly)
good_roots = np.real([r for r in roots if (np.abs(r) == r) and (0 < r < 1)])
np.sort(good_roots)
return good_roots[0]
def get_boundary_gamma(lbda):
lbda = np.atleast_1d(lbda)
dm = np.array([find_boundary_roots(l) for l in lbda])
gamma = -lbda/(2 * mu * dm) + 3 * dm * (1 - dm**3) / mu**2
return gamma
### Compute the boundary line
lbda_max = 6 * (2 / 5.)**(2/3.) * 0.6 / mu
gamma_max = 3 * 4**(-1 / 3.) * 0.75 / mu**2
b_lbdas = np.linspace(0, lbda_max, 100)
b_gammas = get_boundary_gamma(b_lbdas)
print('''Maximum value for Lambda: %.3f \n'''
'''Maximum value for Gamma: %.3f''' % (lbda_max, gamma_max))
Maximum value for Lambda: 0.525 Maximum value for Gamma: 0.102
### Compute the value of delta over a grid in the (Lambda, Gamma) plane
lbdas, gammas = np.meshgrid(np.linspace(0, lbda_max * 1.1, 512),
np.linspace(0, gamma_max * 1.5, 512))
roots = []
for l, g in zip(lbdas.ravel(), gammas.ravel()):
roots.append(find_grad_roots(l, g))
roots = np.array(roots).reshape(gammas.shape)
### 2D case
gammas_2D = np.linspace(0, 2/mu**2, 20)
lambdas_max_2D = ((4 - 2 * gammas_2D * mu**2) / 3.)**(3./2.) / mu
### Plot everything
fig, ax = plt.subplots()
roots[np.isnan(roots)] = 0.5
lbdas, gammas = np.meshgrid(np.linspace(0, lbda_max * 1.1, 512), np.linspace(0, gamma_max * 1.5, 512))
contour_set = ax.contourf(lbdas, gammas, roots.clip(0.5, 1), 256, cmap='gray')
fig.colorbar(contour_set, ticks=[0.5, 0.6, 0.7, 0.8, 0.9, 1.])
boundary = ax.plot(b_lbdas, b_gammas, 'g-', lw=3, alpha=0.8, label='3D case')
boundary2D = ax.plot(lambdas_max_2D, gammas_2D, 'b-',
lw=3, alpha=0.8, label='2D case')
ax.plot(lbda_case1, gamma_case1, 'ko')
ax.set_xlabel(r'Line tension $\bar\Lambda$')
ax.set_ylabel(r'Contractility $\bar\Gamma$')
ax.set_title('Values of $\delta$')
ax.legend(loc='upper right')
fig.savefig(lj.data.get_image('phase_space_boundaries.svg'))
We can now get the values of the optimal scaling $\delta_o$ and of the ground state energy $E_0$
with a simple root finding on the isotropic_grad
function:
delta_o = find_grad_roots(lbda_case1, gamma_case1)
print("Optimum scaling factor = %.3f" % delta_o)
ground_energy = isotropic_energy(delta_o, gamma_case1, lbda_case1)
print("Ground energy = %.3f" % ground_energy)
Optimum scaling factor = 0.887 Ground energy = 0.462
Let's note $A_i$ and $\rho_i$ the initial area and radius, before relaxation, and $A_f, \rho_f$ the same variables after the isotropic relaxation is applied. We now want to multiply all the distances by a factor $c$ such that $h_fA_f = \delta_o^3h_0A_0$, with $\rho_f = c\rho_i$ and $A_f = c^2A_i$.
By noting $\delta_i = h_i / h_0$, we simply have $c = \delta_o / \delta_i$
We now have the two phase space boundaries given by Faradhifar on figure 1. When we pass to the simulation we'll use all that to ensure that optimisation will converge and that we start the simultations close to the ground state.
In order to find the local energy minimum, as shown below, we use a conjugate gradient method
(namely the fmin_lbfgs
algorithm from SciPy). So we need to compute said gradient at each junction vertex that concur to the optimization.
In order to compute the local minimum of the energy, we calculate its gradient at junction vertex $i$ whose position in the 3D space is defined by the vector $\mathbf{r_i}$ pointing from the coordinate system's origin. An edge between two vertices is noted $\mathbf{r_{ij}} = \mathbf{r_j} - \mathbf{r_i}$. The coefficients noted $c_{ij}$ correspond to the terms of the graph's adjacency matrix, such that they are equal to 1 if an edge exists between vertices $i$ and $j$, and 0 otherwise.
$$ \begin{aligned} \mathbf{\nabla_i} E &= (\frac{\partial E}{\partial x}, \frac{\partial E}{\partial y}, \frac{\partial E}{\partial z}) \\ \mathbf{\nabla_i} E &= \sum_\alpha \left(K (V_\alpha - V_0) \mathbf{\nabla_i} V_\alpha + \Gamma L_\alpha \mathbf{\nabla_i} L_\alpha \right)c_{i \alpha} + \sum_i \Lambda_{ij} \mathbf{\nabla_i} \ell_{ij}c_{ij} \end{aligned} $$We have $$ \begin{aligned} \mathbf{\nabla_i}\ell_{ij} &= \frac{\mathbf{r}{ij}}{\ell{ij}}c_{ij}\ \mathbf{\nabla_i}L_\alpha &= \sum_{kn} \mathbf{\nabla_i} \ell_{kn} c_{\alpha k} c_{\alpha n} = \sum_{j} \mathbf{\nabla_i}\ell_{ij} c_{ij} c_{\alpha i} c_{\alpha j} = \sum_{j} \frac{\mathbf{r_{ij}}}{\ell_{ij}}c_{ij}c_{\alpha i} c_{\alpha j}\ \mathbf{\nabla_i}V_{\alpha} &= \sum_{km}\mathbf{\nabla_i}V_{\alpha km} c_{\alpha k}c_{\alpha m}c_{km} = \sum_{km} \mathbf{\nabla_i}h_\alpha A_{\alpha km} c_{\alpha k}c_{\alpha m}c_{km} = A_\alpha \mathbf{\nabla_i}h_\alpha + h_\alpha \sum_{km} \mathbf{\nabla_i} A_{\alpha km}c_{\alpha k} c_{\alpha m}c_{km}\ \mathbf{\nabla_i}V_{\alpha} &= \frac{h_\alpha}{2} \sum_{km}\mathbf{u}{\alpha km} \times \frac{\mathbf{r}{km}}{\nu}c_{\alpha k}c_{\alpha m}c_{km}
\end{aligned} $$
With the normal perpendicular to the $(\mathbf{r_k}, \mathbf{r_m})$ plane noted:
$$ \mathbf{u}_{\alpha km} = \frac{\mathbf{r}_{\alpha k} \times \mathbf{r}_{\alpha m}}{||\mathbf{r}_{\alpha k} \times \mathbf{r}_{\alpha m}||} $$!ls
0 - Introduction.ipynb 5 - Joint.ipynb nb_init.py 1 - Epithelium architecture.ipynb 8 - Parameters exploration.ipynb sandbox.ipynb 2 - Energy minimisation.ipynb Analysis.ipynb saved_graphs 3 - Cell division and type 1 transition.ipynb Figures.ipynb 4 - Single cell apoptosis.ipynb latest.ipynb
import tempfile
d = tempfile.mkdtemp()
d
'/tmp/tmpgyfdubak'
### The script `nb_init` does all the imports
%run nb_init.py
## Epithelium instantiation
eptm = lj.Epithelium(graphXMLfile=lj.data.small_xml(),
save_dir='saved_graphs',
identifier='tuto',
copy=True)
## Scaling so that we start at global epithelium
#eptm.isotropic_relax()
2015-01-26 10:07:01,023 -leg_joint -INFO -successfully imported leg_joint 2015-01-26 10:07:01,024 -leg_joint.epithelium -INFO -Instanciating epithelium tuto 2015-01-26 10:07:01,075 -leg_joint.epithelium -INFO -Initial cells 2015-01-26 10:07:01,075 -leg_joint.epithelium -INFO -Initial junctions 2015-01-26 10:07:01,319 -leg_joint.epithelium -INFO -Update geometry
To do so, we use graph-tool's filtering mechanism. We use two attributes of the Epithelium
object:
is_cell_vert
and is_alive
, which are graph-tools boolean PropertyMaps
.
## Create a filter with only the 'live' cell vertices
vfilt = eptm.is_cell_vert.copy()
vfilt.a *= eptm.is_alive.a
## Filtering
eptm.graph.set_vertex_filter(vfilt)
area = eptm.cells.areas.fa.mean()
height = eptm.rhos.fa.mean() - eptm.rho_lumen
eptm.delta_o = eptm.find_grad_roots()
print('Square root of the ratio between the prefered area\n'
'and the actual average area: %.6f\n'
% (np.sqrt(area / eptm.params["prefered_area"])))
print('Ratio between the prefered height\n'
'and the actual average height: %.6f\n'
% (height / eptm.params['prefered_height']))
print('Scaling factor with respect to elastic volume equilibrium: %.6f\n'
% eptm.delta_o)
# Remove the filter
eptm.graph.set_vertex_filter(None)
Square root of the ratio between the prefered area and the actual average area: 0.886593 Ratio between the prefered height and the actual average height: 0.886593 Scaling factor with respect to elastic volume equilibrium: 0.886593
All three numbers should be equal by definition
The graph bellow corresponds to the analytical computations above
fig, ax = plt.subplots(figsize=(6, 4))
lbda = eptm.paramtree.relative_dic['line_tension']
gamma = eptm.paramtree.relative_dic['contractility']
deltas = np.linspace(0.4, 1.2)
ax.plot(deltas, eptm.isotropic_energy(deltas), 'k-',
label='Analytical total')
ax.plot(eptm.delta_o, eptm.isotropic_energy(eptm.delta_o), 'ro')
ax.set_xlabel('Normalized energy')
ax.set_ylabel('Scaling');
deltas = np.linspace(0.5, 1.2, 20)
fig, ax = plt.subplots(figsize=(5, 8))
lbda = eptm.paramtree.relative_dic['line_tension']
gamma = eptm.paramtree.relative_dic['contractility']
ax.plot(deltas, eptm.isotropic_energy(deltas), 'k-',
label='Analytical total')
ax.plot(eptm.delta_o, eptm.isotropic_energy(eptm.delta_o), 'ro')
ax.plot(deltas, elasticity(deltas), 'b-',
label='Analytical volume elasticity')
ax.plot(deltas, contractility(deltas, gamma), c='orange', ls='-',
label='Analytical contractility')
ax.plot(deltas, tension(deltas, lbda), 'g-',
label='Analytical line tension')
ax.set_xlabel(r'Isotropic scaling $\delta$')
ax.set_ylabel(r'Isotropic energie $\bar E$')
eptm.update_rhotheta()
area0 = eptm.params['prefered_area']
h_0 = eptm.params['prefered_height']
### Cells only area and height
eptm.set_vertex_state([(eptm.is_cell_vert, False),
(eptm.is_alive, False)])
area_avg = eptm.cells.areas.fa.mean()
rho_avg = eptm.rhos.fa.mean()
eptm.set_vertex_state()
### Set height and area to height0 and area0
scale = (area0 / area_avg)**0.5
eptm.scale(scale)
eptm.rho_lumen = rho_avg * scale - h_0
eptm.update_geometry()
energy = eptm.calc_energy()
norm = (eptm.params['prefered_area']**2
* eptm.params['prefered_height']**2
* eptm.params['vol_elasticity']
* (eptm.is_cell_vert.a * eptm.is_alive.a).sum())
cell = eptm.graph.vertex(60)
energies = np.zeros((deltas.size, 3))
#scales = np.linspace(0.5, 1.2, 20) / eptm.delta_o
for n, delta in enumerate(deltas):
eptm.scale(delta)
eptm.update_geometry()
Ec, Ev = eptm.calc_cells_energy()
El, Er = eptm.calc_junctions_energy()
energies[n, :] = [Ev.sum(), El.sum(), Ec.sum()]
eptm.scale(1 / delta)
eptm.isotropic_relax()
energies = np.array(energies) / norm
ax.plot(deltas, energies[:, 0], 'bo:', lw=2, alpha=0.8,
label='Computed volume elasticity')
ax.plot(deltas, energies[:, 1], 'go:', lw=2, alpha=0.8,
label='Computed line tension')
ax.plot(deltas, energies[:, 2], ls=':',
marker='o', c='orange', label='Computed contractility')
ax.plot(deltas, energies.sum(axis=1), 'ko:', lw=2, alpha=0.8,
label='Computed total')
ax.legend(loc='upper left', fontsize=10)
ax.set_ylim(0, 1.2)
plt.savefig(lj.data.get_image('scaling_comparison.svg'))
print(eptm.delta_o, deltas[energies.sum(axis=1).argmin()] * eptm.delta_o)
0.886592687387 0.769935754836
Note that at initialisation time, we also
invoked the isotropic_relax
method,
and we start with an hexagonal lattice, so we should be fairly close to
the ground state energy, except from the cells closer to edges along the $z$ axis.