# Quasistatic vertex models¶

Quasitatic vertex models are based on the minimisation of a potential energy depending on the tissue geometry.

It is the original formalism defined by Honda et al., and in its most widely used form by Farhadifar et al..

tyssue provides a way to easily define such a potential and find the corresponding energy minimum through gradient descent.

:

import sys
import pandas as pd
import numpy as np
import json
import matplotlib.pylab as plt
%matplotlib inline

from scipy import optimize

from tyssue import Sheet, config
from tyssue import SheetGeometry as geom
from tyssue.dynamics import effectors, model_factory
from tyssue.dynamics import SheetModel as model
from tyssue.solvers.quasistatic import QSSolver
from tyssue.draw import sheet_view
from tyssue.draw.plt_draw import plot_forces
from tyssue.io import hdf5


:

h5store = 'data/small_hexagonal.hf5'
#h5store = 'data/before_apoptosis.hf5'

data_names=['face', 'vert', 'edge'])
specs = config.geometry.cylindrical_sheet()
sheet = Sheet('emin', datasets, specs)
sheet.sanitize(trim_borders=True, order_edges=True)

geom.update_all(sheet)

fig, ax = sheet_view(sheet, ['z', 'x'], mode='quick') ## Effector classes¶

An effector designates a dynamical term in the epithelium governing equation. For quasistatic models, we need to provide a method to compute the energy associated with this effector, and the corresponding gradient.

For example, we can consider a line tension effector. The energy is $$E_t = \sum_{ij} \Lambda\ell_{ij}$$ where the sum is over all the edges. For each half-edge, the gradient of $$E_t$$ is defined by two terms, one for the gradient term associated with the half-edge $${ij}$$ source, the over for it’s target. For the $$x$$ coordinate:

$\frac{\partial E_t}{\partial x_i} = \Lambda\left(\sum_k \frac{\partial \ell_{ik}}{\partial x_i} + \sum_m \frac{\partial \ell_{mi}}{\partial x_i}\right) = \Lambda\left(\sum_k \frac{x_i}{\ell_{ik}} - \sum_m \frac{x_i}{\ell_{mi}}\right)$

Where $$\sum_k$$ is are over all the edges which vertex $$i$$ is a source, and $$\sum_m$$ over all the edges of which vertex i is a target.

Here is the definition of the line tension effector:

class LineTension(AbstractEffector):

dimensions = units.line_tension
magnitude = "line_tension"
label = "Line tension"
element = "edge"
specs = {"edge": {"is_active": 1, "line_tension": 1e-2}}

spatial_ref = "mean_length", units.length

@staticmethod
def energy(eptm):
return eptm.edge_df.eval(
"line_tension" "* is_active" "* length / 2"
)  # accounts for half edges

@staticmethod
eptm.edge_df.eval("line_tension * is_active/2"), len(eptm.coords)
)
grad_srce.columns = ["g" + u for u in eptm.coords]


For quasistatic methods, this effectors is agregated with others to define a model object. This object will have two methods compute_energy and compute_gradient that take an Epithelium object as single argument.

:

model_specs = {
"edge": {
"line_tension": 1.,
},
"face": {
"contractility": 1.,
"vol_elasticity": 1.,
"prefered_vol": sheet.face_df["vol"].mean(),
},
}

sheet.update_specs(model_specs, reset=True)

Reseting column line_tension of the edge dataset with new specs
Reseting column contractility of the face dataset with new specs
Reseting column vol_elasticity of the face dataset with new specs
Reseting column prefered_vol of the face dataset with new specs

:

Et, Ec, Ev = model.compute_energy(sheet, full_output=True)
energy = model.compute_energy(sheet, full_output=False)
print('Total energy: {}'.format(energy))

Total energy: 15221.226320514668

:

grad_E = model.compute_gradient(sheet)


/home/docs/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.7/site-packages/tyssue/utils/utils.py:35: FutureWarning: Support for multi-dimensional indexing (e.g. obj[:, None]) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
df_nd = df[:, None]#np.asarray(df).repeat(ndim).reshape((df.size, ndim))

:

gx gy gz
srce
0 -377.344701 -1.563013e+02 697.693957
1 -377.344701 1.563013e+02 697.693957
2 -259.407529 1.074501e+02 -622.859136
3 125.420612 -3.979039e-13 -617.293220
4 -259.407529 -1.074501e+02 -622.859136
:

model.labels

:

['Line tension', 'Contractility', 'Volume elasticity']

:

gradients = model.compute_gradient(sheet, components=True)
gradients = {label: (srce, trgt) for label, (srce, trgt)

/home/docs/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.7/site-packages/tyssue/utils/utils.py:35: FutureWarning: Support for multi-dimensional indexing (e.g. obj[:, None]) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
df_nd = df[:, None]#np.asarray(df).repeat(ndim).reshape((df.size, ndim))

:

gx gy gz
edge
0 -0.00000 -0.000000 -0.500000
1 -0.00000 -0.000000 0.500000
2 0.08434 0.424007 0.251207
3 -0.08434 0.424007 -0.251207
4 0.08434 -0.424007 0.251207
:

%pdb

Automatic pdb calling has been turned ON

:

solver = QSSolver()

scale = 0.2
fig, ax = plot_forces(sheet, geom, model, ['z', 'x'], scale)
fig.set_size_inches(10, 12)
for n, (vx, vy, vz) in sheet.vert_df[sheet.coords].iterrows():
shift = 0.6 * np.sign(vy)
ax.text(vz+shift-0.3, vx, str(n))

def draw_approx_force(ax=None):
fig, ax = plot_forces(sheet, geom, model,
['z', 'x'], scaling=scale, ax=ax,
fig.set_size_inches(10, 12)
return fig, ax

## Uncomment bellow to recompute
fig, ax = draw_approx_force(ax=ax)
#fig

/home/docs/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.7/site-packages/tyssue/utils/utils.py:35: FutureWarning: Support for multi-dimensional indexing (e.g. obj[:, None]) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
df_nd = df[:, None]#np.asarray(df).repeat(ndim).reshape((df.size, ndim)) :




/home/docs/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.7/site-packages/tyssue/utils/utils.py:35: FutureWarning: Support for multi-dimensional indexing (e.g. obj[:, None]) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
df_nd = df[:, None]#np.asarray(df).repeat(ndim).reshape((df.size, ndim))

Error on the gradient (non-dim, per vertex): 1.970e-02

:

options ={'disp':False,
'ftol':1e-5,
'gtol':1e-5}

res = solver.find_energy_min(sheet, geom, model, options=options)
print(res['success'])

True

:

res['message']

:

b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'

:

res['fun']/sheet.face_df.is_alive.sum()

:

240.67238964483766

:

fig, ax = plot_forces(sheet, geom, model, ['z', 'y'], 1)
fig.set_size_inches(10, 12) [ ]: