Packing and influence on cell motility#
Some time ago Hadrien Mary @Hadim_ asked me if tyssue could reproduce this work:
The subtle mechanics of densely packed cells may help explain why some cancerous tumors stay put while others spread https://t.co/G6lbhLhQxQ pic.twitter.com/AK9SzzUazU
— Quanta Magazine (@QuantaMagazine) 11 septembre 2017
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
/tmp/ipykernel_3713/1783676346.py:2: DeprecationWarning:
Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),
(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)
but was not found to be installed on your system.
If this would cause problems for you,
please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466
import pandas as pd
collision solver could not be imported You may need to install CGAL and re-install tyssue
This module needs ipyvolume to work.
You can install it with:
$ conda install -c conda-forge ipyvolume
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")
CGAL-based mesh generation utilities not found, you may need to install CGAL and build from source
C++ extensions are not available for this version
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)
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)
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
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)
/tmp/ipykernel_3713/1283920589.py in ?()
----> 1 solver.solve(tf=120.0, dt=0.1)
2
3 highlight_faces(sheet.face_df, [pulled,], reset_visible=True)
4
~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/tyssue/solvers/viscous.py in ?(self, tf, dt, on_topo_change, topo_change_args)
120 pos = pos + dot_r * dt
121 self.set_pos(pos)
122 self.prev_t = t
123 if self.manager is not None:
--> 124 self.manager.execute(self.eptm)
125 self.geom.update_all(self.eptm)
126 self.manager.update()
127
~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/tyssue/behaviors/event_manager.py in ?(self, eptm)
109 elem_id = kwargs["elem_id"]
110 else:
111 elem_id = -1
112 logger.debug(f"{self.clock}, {elem_id}, {behavior.__name__}")
--> 113 behavior(eptm, self, **kwargs)
~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/tyssue/utils/decorators.py in ?(*args, **kwargs)
62 face = sheet.idx_lookup(face_id, "face")
63 if face is None:
64 return
65 kwargs["face"] = face
---> 66 return func(*args, **kwargs)
/tmp/ipykernel_3713/3300601360.py in ?(sheet, manager, **kwargs)
40 def traction(sheet, manager, **kwargs):
41
42 traction_spec = default_traction_spec
43 traction_spec.update(**kwargs)
---> 44 pulling = tract(sheet, **traction_spec)
45
46 manager.append(traction, **traction_spec)
/tmp/ipykernel_3713/3300601360.py in ?(eptm, face, pull_axis, value, pull_column, face_id)
18 v_edges = eptm.edge_df.query(
19 f'(srce == {pulling_vert}) & (face != {face})'
20 )
21 pulling_edges = v_edges[~v_edges['trgt'].isin(verts)].index
---> 22 eptm.edge_df[pull_column] = 0
23
24 if pulling_edges.size:
25 eptm.edge_df.loc[pulling_edges, pull_column] = value
~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/pandas/core/frame.py in ?(self, key, value)
4295 # Column to set is duplicated
4296 self._setitem_array([key], value)
4297 else:
4298 # set column
-> 4299 self._set_item(key, value)
~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/pandas/core/frame.py in ?(self, key, value)
4522 if isinstance(existing_piece, DataFrame):
4523 value = np.tile(value, (len(existing_piece.columns), 1)).T
4524 refs = None
4525
-> 4526 self._set_item_mgr(key, value, refs)
~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/pandas/core/frame.py in ?(self, key, value, refs)
4474 except KeyError:
4475 # This item wasn't present, just insert at end
4476 self._mgr.insert(len(self._info_axis), key, value, refs)
4477 else:
-> 4478 self._iset_item_mgr(loc, value, refs=refs)
4479
4480 # check if we are modifying a copy
4481 # try to set first as we want an invalid
~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/pandas/core/frame.py in ?(self, loc, value, inplace, refs)
4463 refs: BlockValuesRefs | None = None,
4464 ) -> None:
4465 # when called from _set_item_mgr loc can be anything returned from get_loc
4466 self._mgr.iset(loc, value, inplace=inplace, refs=refs)
-> 4467 self._clear_item_cache()
~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/pandas/core/frame.py in ?(self)
-> 4611 def _clear_item_cache(self) -> None:
4612 self._item_cache.clear()
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")
