"""
Generic forces and energies
"""
import numpy as np
from ..utils import to_nd
from . import units
from .bulk_gradients import lumen_volume_grad, volume_grad
from .planar_gradients import area_grad as area_grad2d
from .planar_gradients import lumen_area_grad
from .sheet_gradients import area_grad, height_grad
[docs]def elastic_force(element_df, var, elasticity, prefered):
params = {"x": var, "K": elasticity, "x0": prefered}
force = element_df.eval("{K} * ({x} - {x0})".format(**params))
return force
def _elastic_force(element_df, x, elasticity, prefered):
force = element_df[elasticity] * (element_df[x] - element_df[prefered])
return force
[docs]def elastic_energy(element_df, var, elasticity, prefered):
params = {"x": var, "K": elasticity, "x0": prefered}
energy = element_df.eval("0.5 * {K} * ({x} - {x0}) ** 2".format(**params))
return energy
def _elastic_energy(element_df, x, elasticity, prefered):
energy = 0.5 * element_df[elasticity] * (element_df[x] - element_df[prefered]) ** 2
return np.array(energy)
[docs]class AbstractEffector:
"""The effector class is used by model factories
to construct a model.
"""
dimensions = None
magnitude = None
spatial_ref = None, None
temporal_ref = None, None
label = "Abstract effector"
element = None # cell, face, edge or vert
specs = {"cell": {}, "face": {}, "edge": {}, "vert": {}}
[docs] @staticmethod
def energy(eptm):
raise NotImplementedError
[docs] @staticmethod
def gradient(eptm):
raise NotImplementedError
[docs] @staticmethod
def get_nrj_norm(specs):
raise NotImplementedError
# @classmethod
# @property
# def __doc__(cls):
# f"""Effector implementing {cls.label} with a magnitude factor
# {cls.magnitude}.
# Works on an `Epithelium` object's {cls.element} elements.
# """
[docs]class LengthElasticity(AbstractEffector):
"""Elastic half edge"""
dimensions = units.line_elasticity
label = "Length elasticity"
magnitude = "length_elasticity"
element = "edge"
spatial_ref = "prefered_length", units.length
specs = {
"edge": {
"is_active": 1,
"length": 1.0,
"length_elasticity": 1.0,
"prefered_length": 1.0,
"ux": (1 / 3) ** 0.5,
"uy": (1 / 3) ** 0.5,
"uz": (1 / 3) ** 0.5,
}
}
[docs] @staticmethod
def get_nrj_norm(specs):
return (
specs["edge"]["length_elasticity"] * specs["edge"]["prefered_length"] ** 2
)
[docs] @staticmethod
def energy(eptm):
return elastic_energy(
eptm.edge_df, "length", "length_elasticity * is_active", "prefered_length"
)
[docs] @staticmethod
def gradient(eptm):
kl_l0 = elastic_force(
eptm.edge_df, "length", "length_elasticity * is_active", "prefered_length"
)
grad = eptm.edge_df[eptm.ucoords] * to_nd(kl_l0, eptm.dim)
grad.columns = ["g" + u for u in eptm.coords]
return -grad, grad
[docs]class PerimeterElasticity(AbstractEffector):
"""From Mapeng Bi et al. https://doi.org/10.1038/nphys3471"""
dimensions = units.line_elasticity
magnitude = "perimeter_elasticity"
label = "Perimeter Elasticity"
element = "face"
specs = {
"face": {
"is_alive": 1,
"perimeter": 1.0,
"perimeter_elasticity": 0.1,
"prefered_perimeter": 3.81,
}
}
spatial_ref = "prefered_perimeter", units.length
[docs] @staticmethod
def energy(eptm):
return eptm.face_df.eval(
"0.5 * is_alive"
"* perimeter_elasticity"
"* (perimeter - prefered_perimeter)** 2"
)
[docs] @staticmethod
def gradient(eptm):
gamma_ = eptm.face_df.eval(
"perimeter_elasticity * is_alive" "* (perimeter - prefered_perimeter)"
)
gamma = eptm.upcast_face(gamma_)
grad_srce = -eptm.edge_df[eptm.ucoords] * to_nd(gamma, len(eptm.coords))
grad_srce.columns = ["g" + u for u in eptm.coords]
grad_trgt = -grad_srce
return grad_srce, grad_trgt
[docs]class FaceAreaElasticity(AbstractEffector):
dimensionless = False
dimensions = units.area_elasticity
magnitude = "area_elasticity"
label = "Area elasticity"
element = "face"
specs = {
"face": {
"is_alive": 1,
"area": 1.0,
"area_elasticity": 1.0,
"prefered_area": 1.0,
},
"edge": {"sub_area": 1 / 6.0},
}
spatial_ref = "prefered_area", units.area
[docs] @staticmethod
def get_nrj_norm(specs):
return specs["face"]["area_elasticity"] * specs["face"]["prefered_area"] ** 2
[docs] @staticmethod
def energy(eptm):
return elastic_energy(
eptm.face_df, "area", "area_elasticity * is_alive", "prefered_area"
)
[docs] @staticmethod
def gradient(eptm):
ka_a0_ = elastic_force(
eptm.face_df, "area", "area_elasticity * is_alive", "prefered_area"
)
ka_a0 = to_nd(eptm.upcast_face(ka_a0_), len(eptm.coords))
if len(eptm.coords) == 2:
grad_a_srce, grad_a_trgt = area_grad2d(eptm)
elif len(eptm.coords) == 3:
grad_a_srce, grad_a_trgt = area_grad(eptm)
grad_a_srce = ka_a0 * grad_a_srce
grad_a_trgt = ka_a0 * grad_a_trgt
grad_a_srce.columns = ["g" + u for u in eptm.coords]
grad_a_trgt.columns = ["g" + u for u in eptm.coords]
return grad_a_srce, grad_a_trgt
[docs]class FaceVolumeElasticity(AbstractEffector):
dimensions = units.vol_elasticity
magnitude = "vol_elasticity"
label = "Volume elasticity"
element = "face"
specs = {
"face": {"is_alive": 1, "vol": 1.0, "vol_elasticity": 1.0, "prefered_vol": 1.0},
"vert": {"height": 1.0},
"edge": {"sub_area": 1 / 6},
}
spatial_ref = "prefered_vol", units.vol
[docs] @staticmethod
def get_nrj_norm(specs):
return specs["face"]["vol_elasticity"] * specs["face"]["prefered_vol"] ** 2
[docs] @staticmethod
def energy(eptm):
return elastic_energy(
eptm.face_df, "vol", "vol_elasticity * is_alive", "prefered_vol"
)
[docs] @staticmethod
def gradient(eptm):
kv_v0_ = elastic_force(
eptm.face_df, "vol", "vol_elasticity * is_alive", "prefered_vol"
)
kv_v0 = to_nd(eptm.upcast_face(kv_v0_), 3)
edge_h = to_nd(eptm.upcast_srce(eptm.vert_df["height"]), 3)
area_ = eptm.edge_df["sub_area"]
area = to_nd(area_, 3)
grad_a_srce, grad_a_trgt = area_grad(eptm)
grad_h = eptm.upcast_srce(height_grad(eptm))
grad_v_srce = kv_v0 * (edge_h * grad_a_srce + area * grad_h)
grad_v_trgt = kv_v0 * (edge_h * grad_a_trgt)
grad_v_srce.columns = ["g" + u for u in eptm.coords]
grad_v_trgt.columns = ["g" + u for u in eptm.coords]
return grad_v_srce, grad_v_trgt
[docs]class CellAreaElasticity(AbstractEffector):
dimensions = units.area_elasticity
magnitude = "area_elasticity"
label = "Area elasticity"
element = "cell"
specs = {
"cell": {
"is_alive": 1,
"area": 1.0,
"area_elasticity": 1.0,
"prefered_area": 1.0,
}
}
spatial_ref = "prefered_area", units.area
[docs] @staticmethod
def get_nrj_norm(specs):
return specs["cell"]["area_elasticity"] * specs["cell"]["prefered_area"] ** 2
[docs] @staticmethod
def energy(eptm):
return elastic_energy(eptm.cell_df, "area", "area_elasticity", "prefered_area")
[docs] @staticmethod
def gradient(eptm):
ka_a0_ = elastic_force(
eptm.cell_df, "area", "area_elasticity * is_alive", "prefered_area"
)
ka_a0 = to_nd(eptm.upcast_cell(ka_a0_), 3)
grad_a_srce, grad_a_trgt = area_grad(eptm)
grad_a_srce = ka_a0 * grad_a_srce
grad_a_trgt = ka_a0 * grad_a_trgt
grad_a_srce.columns = ["g" + u for u in eptm.coords]
grad_a_trgt.columns = ["g" + u for u in eptm.coords]
return grad_a_srce, grad_a_trgt
[docs]class CellVolumeElasticity(AbstractEffector):
dimensions = units.vol_elasticity
magnitude = "vol_elasticity"
label = "Volume elasticity"
element = "cell"
spatial_ref = "prefered_vol", units.vol
specs = {
"cell": {"is_alive": 1, "vol": 1.0, "vol_elasticity": 1.0, "prefered_vol": 1.0}
}
[docs] @staticmethod
def get_nrj_norm(specs):
return specs["cell"]["vol_elasticity"] * specs["cell"]["prefered_vol"] ** 2
[docs] @staticmethod
def energy(eptm):
return elastic_energy(eptm.cell_df, "vol", "vol_elasticity", "prefered_vol")
[docs] @staticmethod
def gradient(eptm):
kv_v0_ = elastic_force(
eptm.cell_df, "vol", "vol_elasticity * is_alive", "prefered_vol"
)
kv_v0 = to_nd(eptm.upcast_cell(kv_v0_), 3)
grad_v_srce, grad_v_trgt = volume_grad(eptm)
grad_v_srce = kv_v0 * grad_v_srce
grad_v_trgt = kv_v0 * grad_v_trgt
grad_v_srce.columns = ["g" + u for u in eptm.coords]
grad_v_trgt.columns = ["g" + u for u in eptm.coords]
return grad_v_srce, grad_v_trgt
[docs]class LumenVolumeElasticity(AbstractEffector):
"""
Global volume elasticity of the object.
For example the volume of the yolk in the Drosophila embryo
"""
dimensions = units.vol_elasticity
magnitude = "lumen_vol_elasticity"
label = "Lumen volume elasticity"
element = "settings"
spatial_ref = "lumen_prefered_vol", units.vol
specs = {
"settings": {
"lumen_vol": 1.0,
"lumen_vol_elasticity": 1.0,
"lumen_prefered_vol": 1.0,
}
}
[docs] @staticmethod
def get_nrj_norm(specs):
return (
specs["settings"]["lumen_vol_elasticity"]
* specs["settings"]["lumen_prefered_vol"] ** 2
)
[docs] @staticmethod
def energy(eptm):
return _elastic_energy(
eptm.settings, "lumen_vol", "lumen_vol_elasticity", "lumen_prefered_vol"
)
[docs] @staticmethod
def gradient(eptm):
kv_v0 = _elastic_force(
eptm.settings, "lumen_vol", "lumen_vol_elasticity", "lumen_prefered_vol"
)
grad_v_srce, grad_v_trgt = lumen_volume_grad(eptm)
grad_v_srce = kv_v0 * grad_v_srce
grad_v_trgt = kv_v0 * grad_v_trgt
grad_v_srce.columns = ["g" + u for u in eptm.coords]
grad_v_trgt.columns = ["g" + u for u in eptm.coords]
return grad_v_srce, grad_v_trgt
[docs]class LineTension(AbstractEffector):
dimensions = units.line_tension
magnitude = "line_tension"
label = "Line tension"
element = "edge"
specs = {"edge": {"is_active": 1, "line_tension": 1.0}}
spatial_ref = "mean_length", units.length
[docs] @staticmethod
def energy(eptm):
return eptm.edge_df.eval(
"line_tension" "* is_active" "* length / 2"
) # accounts for half edges
[docs] @staticmethod
def gradient(eptm):
grad_srce = -eptm.edge_df[eptm.ucoords] * to_nd(
eptm.edge_df.eval("line_tension * is_active/2"), len(eptm.coords)
)
grad_srce.columns = ["g" + u for u in eptm.coords]
grad_trgt = -grad_srce
return grad_srce, grad_trgt
[docs]class FaceContractility(AbstractEffector):
dimensions = units.line_elasticity
magnitude = "contractility"
label = "Contractility"
element = "face"
specs = {"face": {"is_alive": 1, "perimeter": 1.0, "contractility": 1.0}}
spatial_ref = "mean_perimeter", units.length
[docs] @staticmethod
def energy(eptm):
return eptm.face_df.eval("0.5 * is_alive * contractility * perimeter ** 2")
[docs] @staticmethod
def gradient(eptm):
gamma_ = eptm.face_df.eval("contractility * perimeter * is_alive")
gamma = eptm.upcast_face(gamma_)
grad_srce = -eptm.edge_df[eptm.ucoords] * to_nd(gamma, len(eptm.coords))
grad_srce.columns = ["g" + u for u in eptm.coords]
grad_trgt = -grad_srce
return grad_srce, grad_trgt
[docs]class SurfaceTension(AbstractEffector):
dimensions = units.area_tension
magnitude = "surface_tension"
spatial_ref = "prefered_area", units.area
label = "Surface tension"
element = "face"
specs = {"face": {"is_active": 1, "surface_tension": 1.0, "area": 1.0}}
[docs] @staticmethod
def energy(eptm):
return eptm.face_df.eval("surface_tension * area")
[docs] @staticmethod
def gradient(eptm):
G = to_nd(eptm.upcast_face(eptm.face_df["surface_tension"]), len(eptm.coords))
grad_a_srce, grad_a_trgt = area_grad(eptm)
grad_a_srce = G * grad_a_srce
grad_a_trgt = G * grad_a_trgt
grad_a_srce.columns = ["g" + u for u in eptm.coords]
grad_a_trgt.columns = ["g" + u for u in eptm.coords]
return grad_a_srce, grad_a_trgt
[docs]class LineViscosity(AbstractEffector):
dimensions = units.line_viscosity
magnitude = "edge_viscosity"
label = "Linear viscosity"
element = "edge"
spatial_ref = "mean_length", units.length
temporal_ref = "dt", units.time
specs = {"edge": {"is_active": 1, "edge_viscosity": 1.0}}
[docs] @staticmethod
def gradient(eptm):
grad_srce = eptm.edge_df[["vx", "vy", "vz"]] * to_nd(
eptm.edge_df["edge_viscosity"], len(eptm.coords)
)
grad_srce.columns = ["g" + u for u in eptm.coords]
return grad_srce, None
[docs]class BorderElasticity(AbstractEffector):
dimensions = units.line_elasticity
label = "Border edges elasticity"
magnitude = "border_elasticity"
element = "edge"
spatial_ref = "prefered_length", units.length
specs = {
"edge": {
"is_active": 1,
"length": 1.0,
"border_elasticity": 1.0,
"prefered_length": 1.0,
"is_border": 1.0,
}
}
[docs] @staticmethod
def get_nrj_norm(specs):
return (
specs["edge"]["border_elasticity"] * specs["edge"]["prefered_length"] ** 2
)
[docs] @staticmethod
def energy(eptm):
return elastic_energy(
eptm.edge_df,
"length",
"border_elasticity * is_active * is_border / 2",
"prefered_length",
)
[docs] @staticmethod
def gradient(eptm):
kl_l0 = elastic_force(
eptm.edge_df,
var="length",
elasticity="border_elasticity * is_active * is_border",
prefered="prefered_length",
)
grad = eptm.edge_df[eptm.ucoords] * to_nd(kl_l0, eptm.dim)
grad.columns = ["g" + u for u in eptm.coords]
return grad / 2, -grad / 2
[docs]class LumenAreaElasticity(AbstractEffector):
"""
..math: \frac{K_Y}{2}(A_{\mathrm{lumen}} - A_{0,\mathrm{lumen}})^2
"""
dimensions = units.area_elasticity
label = "Lumen volume constraint"
magnitude = "lumen_elasticity"
element = "settings"
spatial_ref = "lumen_prefered_vol", units.area
specs = {
"settings": {
"lumen_elasticity": 1.0,
"lumen_prefered_vol": 1.0,
"lumen_vol": 1.0,
}
}
[docs] @staticmethod
def energy(eptm):
Ky = eptm.settings["lumen_elasticity"]
V0 = eptm.settings["lumen_prefered_vol"]
Vy = eptm.settings["lumen_vol"]
return np.array([Ky * (Vy - V0) ** 2 / 2])
[docs] @staticmethod
def gradient(eptm):
Ky = eptm.settings["lumen_elasticity"]
V0 = eptm.settings["lumen_prefered_vol"]
Vy = eptm.settings["lumen_vol"]
grad_srce, grad_trgt = lumen_area_grad(eptm)
return (Ky * (Vy - V0) * grad_srce, Ky * (Vy - V0) * grad_trgt)
[docs]class RadialTension(AbstractEffector):
"""
Apply a tension perpendicular to a face.
"""
dimensions = units.line_tension
magnitude = "radial_tension"
label = "Apical basal tension"
element = "face"
specs = {"face": {"height": 1.0, "radial_tension": 1.0}}
[docs] @staticmethod
def energy(eptm):
return eptm.face_df.eval("height * radial_tension")
[docs] @staticmethod
def gradient(eptm):
upcast_tension = eptm.upcast_face(
eptm.face_df.eval("radial_tension / num_sides")
)
upcast_height = eptm.upcast_srce(height_grad(eptm))
grad_srce = to_nd(upcast_tension, 3) * upcast_height
grad_srce.columns = ["g" + u for u in eptm.coords]
return grad_srce / 2, grad_srce / 2
[docs]class BarrierElasticity(AbstractEffector):
"""
Barrier use to maintain the tissue integrity.
"""
dimensions = units.line_elasticity
magnitude = "barrier_elasticity"
label = "Barrier elasticity"
element = "vert"
specs = {
"vert": {"barrier_elasticity": 1.0, "is_active": 1, "delta_rho": 0.0}
} # distance to a barrier membrane
[docs] @staticmethod
def energy(eptm):
return eptm.vert_df.eval("delta_rho**2 * barrier_elasticity/2")
[docs] @staticmethod
def gradient(eptm):
grad = height_grad(eptm) * to_nd(
eptm.vert_df.eval("barrier_elasticity * delta_rho"), 3
)
grad.columns = ["g" + c for c in eptm.coords]
return grad, None
[docs]class MidlineBoundary(AbstractEffector):
"""
Elastic boundary at the x-axis, to be used with MidlineBoundaryGeometry.
Intended to use with a high midline_boundary_stiffness to prevent vertices
from crossing the midline.
"""
dimensions = units.line_elasticity
magnitude = "midline_boundary"
label = "Midline boundary"
element = "vert"
specs = {"vert": {"boundary_K": 280, "is_active": 1, "delta_boundary": 0}}
[docs] @staticmethod
def energy(eptm):
return eptm.vert_df.eval(
"0.5 * delta_boundary**2 * {}".format(
str(eptm.settings["midline_boundary_stiffness"])
)
)
[docs] @staticmethod
def gradient(eptm):
# just a bunch of zeros
kl_l0 = elastic_force(eptm.vert_df, "delta_boundary", "0", "0")
grad = eptm.vert_df[eptm.coords] * to_nd(kl_l0, eptm.dim)
grad.columns = ["g" + u for u in eptm.coords]
grad["gy"] = 0
if "z" in eptm.coords:
grad["gz"] = 0
return grad, grad
def _exponants(dimensions, ref_dimensions, spatial_unit=None, temporal_unit=None):
spatial_exponant = time_exponant = 0
rel_dimensionality = (dimensions / ref_dimensions).dimensionality
if spatial_unit is not None:
spatial_exponant = (
rel_dimensionality.get(units.length, 0)
/ spatial_unit.dimensionality[units.length]
)
if temporal_unit is not None:
time_exponant = (
rel_dimensionality.get(units.time, 0)
/ temporal_unit.dimensionality[units.time]
)
return spatial_exponant, time_exponant
[docs]def scaler(nondim_specs, dim_specs, effector, ref_effector):
spatial_val, spatial_unit = ref_effector.spatial_ref
temporal_val, temporal_unit = ref_effector.temporal_ref
s_expo, t_expo = _exponants(
effector.dimensions, ref_effector.dimensions, spatial_unit, temporal_unit
)
ref_magnitude = ref_effector.magnitude
ref_element = ref_effector.element
factor = (
dim_specs[ref_element][ref_magnitude]
* dim_specs[ref_element].get(spatial_val, 1) ** s_expo
* dim_specs[ref_element].get(temporal_val, 1) ** t_expo
)
return factor
[docs]def dimensionalize(nondim_specs, dim_specs, effector, ref_effector):
magnitude = effector.magnitude
element = effector.element
factor = scaler(nondim_specs, dim_specs, effector, ref_effector)
dim_specs[element][magnitude] = factor * nondim_specs[element][magnitude]
return dim_specs
[docs]def normalize(dim_specs, nondim_specs, effector, ref_effector):
magnitude = effector.magnitude
element = effector.element
factor = scaler(nondim_specs, dim_specs, effector, ref_effector)
dim_specs[element][magnitude] = nondim_specs[element][magnitude] / factor
return dim_specs