Reversible Network Reconnection Model#

Attempt at implementing the RNR model as described in Okuda et al. 2012:

Reversible network reconnection model for simulating large deformation in dynamic tissue morphogenesis, Satoru Okuda, Yasuhiro Inoue, Mototsugu Eiraku, Yoshiki Sasai and Taiji Adachi Biomech Model Mechanobiol (2013) 12:627–644 DOI 10.1007/s10237-012-0430-7

The authors give 4 conditions which we detail and implement bellow.

Condition 1 - Center of a face#

The center of a face is defined by the average position of the face’s edges midpoints, weighted by their lengths.

\[ \mathbf{r}_{\alpha} = \frac{\sum_{ij\alpha}\ell_{ij} (\mathbf{r}_i + \mathbf{r}_j)/2}{\sum_{ij\alpha}\ell_{ij}}\]
import numpy as np, pandas as pd
import ipyvolume as ipv
import matplotlib.pyplot as plt
%matplotlib inline

from tyssue import Epithelium
from tyssue import BulkGeometry, RNRGeometry
from tyssue import Sheet
from tyssue.config.geometry import bulk_spec
from tyssue.generation import three_faces_sheet
from tyssue.generation import extrude
from tyssue.draw import sheet_view, highlight_cells

from tyssue.topology.bulk_topology import IH_transition, HI_transition



draw_spec = {'face': {'visible': True}}




sheet = Sheet.planar_sheet_3d('sheet', 5, 5, 1, 1)
sheet.sanitize()
datasets = extrude(sheet.datasets, method='translation')

eptm = Epithelium('20faces_3D', datasets, bulk_spec())
eptm.vert_df[eptm.coords] += np.random.normal(scale=1e-2, size=(eptm.Nv, 3))
RNRGeometry.update_all(eptm)
RNRGeometry.center(eptm)
RNRGeometry.update_all(eptm)


eptm.settings['threshold_length'] = 1e-4
print(eptm.Nf, eptm.Ne, eptm.Nv)
draw_spec['face']['color'] = eptm.face_df.area
draw_spec['face']['color_range'] = (0, 1)

ipv.clear()
fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)
fig
/tmp/ipykernel_3653/2203072209.py:1: 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 numpy as np, pandas as pd
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 2
      1 import numpy as np, pandas as pd
----> 2 import ipyvolume as ipv
      3 import matplotlib.pyplot as plt
      4 get_ipython().run_line_magic('matplotlib', 'inline')

File ~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/ipyvolume/__init__.py:8
      6 from ipyvolume import datasets  # noqa: F401
      7 from ipyvolume import embed  # noqa: F401
----> 8 from ipyvolume.widgets import *  # noqa: F401, F403
      9 from ipyvolume.transferfunction import *  # noqa: F401, F403
     10 from ipyvolume.pylab import *  # noqa: F401, F403

File ~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/ipyvolume/widgets.py:23
     21 import ipyvolume._version
     22 from ipyvolume.traittypes import Image
---> 23 from ipyvolume.serialize import (
     24     array_cube_tile_serialization,
     25     array_serialization,
     26     array_sequence_serialization,
     27     color_serialization,
     28     texture_serialization,
     29 )
     30 from ipyvolume.transferfunction import TransferFunction
     31 from ipyvolume.utils import debounced, grid_slice, reduce_size

File ~/checkouts/readthedocs.org/user_builds/tyssue/conda/latest/lib/python3.12/site-packages/ipyvolume/serialize.py:17
     15 import ipywidgets
     16 import ipywebrtc
---> 17 from ipython_genutils.py3compat import string_types
     19 from ipyvolume import utils
     22 logger = logging.getLogger("ipyvolume")

ModuleNotFoundError: No module named 'ipython_genutils'

I→H transition#

IH transition in a bulk and a RNR

This produces a shape similar to the one studied by LM Escudero and collegues in their article.

e_1011 = 26
eptm.settings['threshold_length'] = 0.2

srce, trgt, face, cell = eptm.edge_df.loc[e_1011,
                                          ['srce', 'trgt',
                                           'face', 'cell']]

IH_transition(eptm, e_1011)
RNRGeometry.update_all(eptm)

highlight_cells(eptm, cell, reset_visible=True)
draw_spec['face']['color'] = eptm.face_df.area


ipv.clear()
fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)
fig
cell 0 is already closed
cell 5 is already closed

H → I transition#

fa = 93


HI_transition(eptm, fa)
RNRGeometry.update_all(eptm)

highlight_cells(eptm, cell, reset_visible=True)
ipv.clear()
draw_spec['face']['color'] = eptm.face_df.area
draw_spec['face']['visible'] = True
eptm.face_df.visible = True
#fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)
#fig
cell 4 is already closed
cell 0 is already closed
cell 5 is already closed
cell 1 is already closed
eptm.edge_df[eptm.edge_df['sub_vol'] <= 1e-10][["srce", "trgt", "face", "cell", "length", "sub_area", "sx", "sy", "sz", "tx", "ty", "tz", "fx", "fy", "fz"]]
srce trgt face cell length sub_area sx sy sz tx ty tz fx fy fz
edge
241 39 51 50 5 0.534805 7.372575e-18 -1.312808 -0.374665 -0.493168 -0.838700 -0.621678 -0.508126 -1.075754 -0.498172 -0.500647
244 51 39 51 0 0.534805 7.564593e-18 -0.838700 -0.621678 -0.508126 -1.312808 -0.374665 -0.493168 -1.075754 -0.498172 -0.500647
fig, ax = sheet_view(eptm, coords=["x", "y"])

ax.plot(eptm.edge_df.loc[241, ["sx", "tx"]], eptm.edge_df.loc[241, ["sy", "ty"]], "ko-")
ax.plot(eptm.edge_df.loc[241, ["fx", "sx"]], eptm.edge_df.loc[241, ["fy", "sy"]], "ro-")

for e, data in eptm.edge_df[eptm.edge_df.face.isin([50, 51])].iterrows():
    
    ax.plot(data[["sx", "tx"]], data[["sy", "ty"]], "k-",lw=0.2)
../_images/6cac387716da851190607d21b70057b92aa3f3ee44b8ab7e7b26c9daed74de3b.png
RNRGeometry.update_areas(eptm)
highlight_cells(eptm, cell, reset_visible=False)
ipv.clear()
draw_spec['face']['color'] = eptm.face_df.area
eptm.face_df.visible = True
#fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)
#fig

Testing for I→H / H →I transition triggers#

highlight_cells(eptm, cell, reset_visible=True)
draw_spec['face']['color'] = eptm.face_df.area


fig, ax = sheet_view(eptm, coords = ["x", "z"], **draw_spec)
fig, ax = sheet_view(eptm, coords = ["x", "y"], **draw_spec)
../_images/77e5aeee653cfef148bd592dd8133c4dd5795492d715f3d1471347c0d6c5f535.png ../_images/c11d210765ab22a8f20f6c6645dfa47221f04aabd454e73915303944e0cb5b9e.png
eptm.settings['threshold_length'] = 1e-3
#eptm.settings['threshold_length'] = eptm.edge_df.length.min()+1e-3
def find_rearangements(eptm):
    l_th = eptm.settings['threshold_length']
    up_num_sides = eptm.upcast_face(eptm.face_df['num_sides'])   
    shorts = eptm.edge_df[eptm.edge_df['length'] < l_th]
    non_triangular = up_num_sides[up_num_sides > 4 ].index
    edges_IH = set(shorts.index).intersection(non_triangular)

    max_f_length = shorts.groupby('face')['length'].apply(max)
    short_faces = max_f_length[max_f_length < l_th].index
    three_faces = eptm.face_df[eptm.face_df['num_sides'] == 3].index
    faces_HI = set(three_faces).intersection(short_faces)
    return edges_IH, faces_HI

find_rearangements(eptm)
(set(), set())

Condition 3#

This condition is satisfied if eptm.settings['threshold_length'] is well defined, i.e, small with respect to the unit length or the average edge length.

Condition 4#

  • (i) Two edges never share two vertices simultaneously.

In our half-edge architecture, this is not as straight forward. But I think the condition can be reformulated as: Two edges from the same face never share two edges simultanously. I don’t know how to demonstrate this, but I think these are equivalent (to get two edges together, you need to “squeeze” one face between those two.

  • (ii) Two polygonal faces never share two or more edges simultaneously.

from tyssue.topology.base_topology import condition_4i, condition_4ii
condition_4i(eptm), condition_4ii(eptm)
(Int64Index([], dtype='int64', name='face'),
 array([], shape=(0, 2), dtype=int64))
eptm.settings['threshold_length'] = 0.5

IH_transition(eptm, eptm.edge_df.index[-1])
RNRGeometry.update_all(eptm)
draw_spec['face']['color'] = eptm.face_df.area

#ipv.clear()
#fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)
#fig
cell 0 is already closed
cell 5 is already closed
/home/guillaume/Dev/tyssue/tyssue/geometry/bulk_geometry.py:107: RuntimeWarning: invalid value encountered in divide
  weighted_pos.values / eptm.face_df["perimeter"].values[:, np.newaxis]
../_images/c11d210765ab22a8f20f6c6645dfa47221f04aabd454e73915303944e0cb5b9e.png