Packing and influence on cell motility

Some time ago Hadrien Mary @Hadim_ asked me if tyssue could reproduce this work:

Cellular traffic jam

Illustration is from Lucy Reading @LucyIkkanda

The relevant model is described in Mapeng Bi et al. (nature physics version ($$))

The master equation is the following:

\[\epsilon = \sum_\alpha \left[ (a_\alpha - 1)^2 + \frac{(p_\alpha - p_0)^2}{r}\right]\]

With a unit prefered area and \(p_\alpha\) the cell perimeter.

Imports

import numpy as np
import pandas as pd

from tyssue import config, Sheet, PlanarGeometry
from tyssue.solvers import QSSolver

from tyssue.solvers.viscous import EulerSolver
from tyssue.draw import sheet_view

from tyssue.dynamics.factory import model_factory
from tyssue.dynamics import effectors, units

from tyssue.utils import to_nd
from tyssue.utils.testing import effector_tester, model_tester
from tyssue.behaviors import EventManager
from tyssue.behaviors.sheet import basic_events

from tyssue.draw import highlight_faces, create_gif

Create a 2D patch of cells

geom = PlanarGeometry

sheet = Sheet.planar_sheet_2d('jam', 15, 15, 1, 1, noise=0.2)
geom.update_all(sheet)

sheet.remove(sheet.cut_out([[0, 8], [0, 8]]), trim_borders=True)
sheet.sanitize()
geom.scale(sheet, sheet.face_df.area.mean()**-0.5, ['x', 'y'])
geom.update_all(sheet)
sheet.reset_index()
sheet.reset_topo()
fig, ax = sheet_view(sheet, mode="quick")
../_images/C-2DMigration_5_0.png

Define the relevant mechanical components

   
# Adding some gigling
class BrownianMotion(effectors.AbstractEffector):
    
    label = 'Brownian Motion'
    element = 'vert'
    specs = {"settings": {"temperature": 1e-3}}
    
    def energy(eptm):
        T = eptm.settings['temperature']
        return np.ones(eptm.Nv) * T / eptm.Nv
    
    def gradient(eptm):
        T = eptm.settings['temperature']
        scale = T/eptm.edge_df.length.mean()
        
        grad = pd.DataFrame(
            data=np.random.normal(0, scale=scale, size=(eptm.Nv, eptm.dim)),
            index=eptm.vert_df.index,
            columns=['g'+c for c in eptm.coords]
        )
        return grad, None
    
    

Quasistatic gradient descent

With only the conservative potential terms

model = model_factory(
    [effectors.PerimeterElasticity,
     effectors.FaceAreaElasticity])


model_specs = {
    'face': {
        'area_elasticity': 1.,
        'prefered_area': 1.,
        'perimeter_elasticity': 0.1, # 1/2r in the master equation
        'prefered_perimeter': 3.81,
        },
    'edge': {
        'ux': 0,
        'uy': 0,
        'is_active': 1,
        'line_tension': 0.0
    },
    'settings': {'temperature': 1e-2}
}
    
sheet.update_specs(model_specs, reset=True)

res = QSSolver().find_energy_min(sheet, PlanarGeometry, model)

print(res.message)
Reseting column ux of the edge dataset with new specs
Reseting column uy of the edge dataset with new specs
CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH

Backup so we can play with parameters

bck = sheet.copy()

Pulling on a face

    
def tract(eptm, face, pull_axis, value, pull_column='line_tension', face_id=None):
    """Modeling face traction as shrinking apical junctions on the neighouring cells,
    
    As if pseudopods where pulling between the cells
    
    
    """
    
    
    pull_axis = np.asarray(pull_axis)
    edges = eptm.edge_df.query(f'face == {face}')
    verts = edges['srce'].values
    r_ai = edges[["r"+c for c in eptm.coords]].values
    proj = (r_ai * pull_axis[None, :]).sum(axis=1)
    
    pulling_vert = verts[np.argmax(proj)]
    
    v_edges = eptm.edge_df.query(
        f'(srce == {pulling_vert}) & (face != {face})'
    )
    pulling_edges = v_edges[~v_edges['trgt'].isin(verts)].index
    eptm.edge_df[pull_column] = 0

    if pulling_edges.size:
        eptm.edge_df.loc[pulling_edges, pull_column] = value
        
        
default_traction_spec = {
    "face": -1,
    "face_id": -1,
    "pull_axis": [0.1, 0.9],
    "value": 4,
    "pull_column": "line_tension"
}


from tyssue.utils.decorators import face_lookup

@face_lookup
def traction(sheet, manager, **kwargs):
    
    traction_spec = default_traction_spec
    traction_spec.update(**kwargs)
    pulling = tract(sheet, **traction_spec)
    
    manager.append(traction, **traction_spec)
pulled = 2

Simple visualisation

highlight_faces(sheet.face_df, [pulled,], reset_visible=True)
fig, ax = sheet_view(
    sheet, 
    mode="2D", 
    face={"visible": True},
    edge={"head_width": 0.0}, 
    vert={"visible": False})
fig.set_size_inches(12, 12)
../_images/C-2DMigration_16_0.png

Model with all the components

model = model_factory(
    [effectors.PerimeterElasticity,
     BrownianMotion,
     effectors.LineTension,
     effectors.FaceAreaElasticity])

Setup the event manager and the solver

sheet = bck.copy()


# setting up values for the whole epithelium
model_specs = {
    'face': {
        'area_elasticity': 1.,
        'prefered_area': 1.,
        'perimeter_elasticity': 0.05, # 1/2r in the master equation
        'prefered_perimeter': 3.81,
        "id": sheet.face_df.index
        },
    'vert': {
        "viscosity": 1.0
        },
    'edge': {'ux': 0, 'uy': 0},
    'settings': {
        'temperature': 2e-1,
        "p_4": 10.0,
        "p_5p": 1.0,
        "threshold_length": 2e-2
    }
}

sheet.update_specs(model_specs, reset=True)

# This allows to auomaticaly solve topology changes

manager = EventManager("face", )
manager.append(basic_events.reconnect)
manager.append(traction, face_id=pulled)

# Implicit Euler solver

solver = EulerSolver(
    sheet,
    geom,
    model,
    manager=manager,
    bounds=(
        -sheet.edge_df.length.median()/10,
        sheet.edge_df.length.median()/10
    )
)
manager.update()


sheet.face_df['prefered_perimeter'] = 3.8
Reseting column area_elasticity of the face dataset with new specs
Reseting column prefered_area of the face dataset with new specs
Reseting column perimeter_elasticity of the face dataset with new specs
Reseting column prefered_perimeter of the face dataset with new specs
Reseting column id of the face dataset with new specs
Reseting column ux of the edge dataset with new specs
Reseting column uy of the edge dataset with new specs

Run the solver

    
solver.solve(tf=120.0, dt=0.1)

highlight_faces(sheet.face_df, [pulled,], reset_visible=True)

fig, ax = sheet_view(
    sheet,
    mode="2D",
    face={"visible": True},
    edge={"head_width": 0.0, "color": sheet.edge_df["line_tension"]},
    vert={"visible": False}
)
fig.set_size_inches(6, 6)
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-11-44dcd6235ffd> in <module>
----> 1 solver.solve(tf=120.0, dt=0.1)
      2 
      3 highlight_faces(sheet.face_df, [pulled,], reset_visible=True)
      4 
      5 fig, ax = sheet_view(

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/tyssue-0.8.1-py3.9-linux-x86_64.egg/tyssue/solvers/viscous.py in solve(self, tf, dt, on_topo_change, topo_change_args)
    125         for t in np.arange(self.prev_t, tf + dt, dt):
    126             pos = self.current_pos
--> 127             dot_r = self.ode_func(t, pos)
    128             if self.bounds is not None:
    129                 dot_r = np.clip(dot_r, *self.bounds)

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/tyssue-0.8.1-py3.9-linux-x86_64.egg/tyssue/solvers/viscous.py in ode_func(self, t, pos)
    156         """
    157 
--> 158         grad_U = self.model.compute_gradient(self.eptm).loc[self.eptm.active_verts]
    159         return (
    160             -grad_U.values

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/tyssue-0.8.1-py3.9-linux-x86_64.egg/tyssue/dynamics/factory.py in compute_gradient(eptm, components)
    116             ]
    117             if trgt_grads:
--> 118                 grad_t = eptm.sum_trgt(sum(trgt_grads))
    119             vert_grads = [g[0] for g in grads if g[0].shape[0] == eptm.Nv]
    120             if vert_grads:

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/tyssue-0.8.1-py3.9-linux-x86_64.egg/tyssue/core/objects.py in sum_trgt(self, df)
    455         summed : :class:`pd.DataFrame` the summed data, indexed by the source vertices.
    456         """
--> 457         return self._lvl_sum(df, "trgt")
    458 
    459     def sum_face(self, df):

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/tyssue-0.8.1-py3.9-linux-x86_64.egg/tyssue/core/objects.py in _lvl_sum(self, df, lvl)
    435             df_ = df.copy()
    436         df_[lvl] = self.edge_df[lvl]
--> 437         return df_.groupby(lvl).sum()
    438 
    439     def sum_srce(self, df):

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/pandas/core/groupby/groupby.py in sum(self, numeric_only, min_count)
   1648         # _agg_general() returns. GH #31422
   1649         with com.temp_setattr(self, "observed", True):
-> 1650             result = self._agg_general(
   1651                 numeric_only=numeric_only,
   1652                 min_count=min_count,

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/pandas/core/groupby/groupby.py in _agg_general(self, numeric_only, min_count, alias, npfunc)
   1022             result = None
   1023             try:
-> 1024                 result = self._cython_agg_general(
   1025                     how=alias,
   1026                     alt=npfunc,

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/pandas/core/groupby/generic.py in _cython_agg_general(self, how, alt, numeric_only, min_count)
   1013         self, how: str, alt=None, numeric_only: bool = True, min_count: int = -1
   1014     ) -> DataFrame:
-> 1015         agg_mgr = self._cython_agg_blocks(
   1016             how, alt=alt, numeric_only=numeric_only, min_count=min_count
   1017         )

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/pandas/core/groupby/generic.py in _cython_agg_blocks(self, how, alt, numeric_only, min_count)
   1116         #  continue and exclude the block
   1117         # NotImplementedError -> "ohlc" with wrong dtype
-> 1118         new_mgr = data.apply(blk_func, ignore_failures=True)
   1119 
   1120         if not len(new_mgr):

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/pandas/core/internals/managers.py in apply(self, f, align_keys, ignore_failures, **kwargs)
    433 
    434         if ignore_failures:
--> 435             return self._combine(result_blocks)
    436 
    437         if len(result_blocks) == 0:

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/pandas/core/internals/managers.py in _combine(self, blocks, copy, index)
    764         for b in blocks:
    765             b = b.copy(deep=copy)
--> 766             b.mgr_locs = inv_indexer[b.mgr_locs.indexer]
    767             new_blocks.append(b)
    768 

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/pandas/core/internals/blocks.py in mgr_locs(self, new_mgr_locs)
    270     def mgr_locs(self, new_mgr_locs):
    271         if not isinstance(new_mgr_locs, libinternals.BlockPlacement):
--> 272             new_mgr_locs = libinternals.BlockPlacement(new_mgr_locs)
    273 
    274         self._mgr_locs = new_mgr_locs

pandas/_libs/internals.pyx in pandas._libs.internals.BlockPlacement.__init__()

~/miniconda3/envs/tyssue/lib/python3.9/site-packages/numpy/core/_asarray.py in require(a, dtype, requirements, like)
    398         requirements.remove('C')
    399 
--> 400     arr = array(a, dtype=dtype, order=order, copy=False, subok=subok)
    401 
    402     for prop in requirements:

KeyboardInterrupt: 

Define a simple drawing function

def view(sheet):
    highlight_faces(sheet.face_df, [pulled], reset_visible=True)
    geom.update_all(sheet)
    if sheet.edge_df['line_tension'].max():
        ecolor = sheet.edge_df['line_tension']
    else:
        ecolor = "#aaaaaa"
        
    fig, ax = sheet_view(
        sheet,
        mode="2D",
        face={"visible": True},
        edge={"head_width": 0.0, "color": ecolor, "width": 2},
        vert={"visible": False}
    )
    fig.set_size_inches(8, 8)
    ax.set_xticks([])
    ax.set_yticks([])
    
    return fig, ax

Create a gif of the resulting simulation

create_gif(solver.history, "demo.gif", draw_func=view, num_frames=120)
solver.history.time_stamps
array([0.000e+00, 1.000e-01, 2.000e-01, ..., 1.198e+02, 1.199e+02,
       1.200e+02])
from IPython.display import Image
Image("demo.gif")
<IPython.core.display.Image object>