A-ReversibleNetworkReconnection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reversible Network Reconnection Model\n",
    "\n",
    "Attempt at implementing the RNR model as described in Okuda et al. 2012:\n",
    "\n",
    ">_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](http://link.springer.com/article/10.1007%2Fs10237-012-0430-7)\n",
    "\n",
    "\n",
    "The authors give 4 conditions which we detail and implement bellow.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Condition 1 - Center of a face\n",
    "\n",
    "The center of a face is defined by the average position of the face's edges midpoints, weighted by their lengths.\n",
    "\n",
    "$$ \\mathbf{r}_{\\alpha} = \\frac{\\sum_{ij\\alpha}\\ell_{ij} (\\mathbf{r}_i + \\mathbf{r}_j)/2}{\\sum_{ij\\alpha}\\ell_{ij}}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "110 468 60\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9207b176f36749e0830e2412127656e9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Figure(camera=PerspectiveCamera(fov=45.0, position=(0.0, 0.0, 2.0), projectionMatrix=(1.0, 0.0, 0.0, 0.0, 0.0,…"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np, pandas as pd\n",
    "import ipyvolume as ipv\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from tyssue import Epithelium\n",
    "from tyssue import BulkGeometry, RNRGeometry, SheetGeometry\n",
    "from tyssue import Sheet\n",
    "from tyssue.utils import to_nd\n",
    "\n",
    "from tyssue.config.geometry import bulk_spec\n",
    "from tyssue.generation import three_faces_sheet\n",
    "from tyssue.generation import extrude\n",
    "from tyssue.draw import sheet_view, highlight_cells\n",
    "\n",
    "from tyssue.topology.bulk_topology import IH_transition, HI_transition\n",
    "from tyssue.topology import find_rearangements\n",
    "\n",
    "from tyssue.dynamics import model_factory, effectors\n",
    "from tyssue.solvers import QSSolver\n",
    "\n",
    "draw_spec = {'face': {'visible': True}}\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "sheet = Sheet.planar_sheet_3d('sheet', 5, 6, 1, 1)\n",
    "sheet.sanitize(trim_borders=True)\n",
    "\n",
    "geom = RNRGeometry()\n",
    "datasets = extrude(sheet.datasets, method='translation')\n",
    "\n",
    "specs = bulk_spec()\n",
    "eptm = Epithelium('20faces_3D', datasets, bulk_spec())\n",
    "RNRGeometry.update_all(eptm)\n",
    "RNRGeometry.center(eptm)\n",
    "RNRGeometry.update_all(eptm)\n",
    "\n",
    "\n",
    "eptm.settings['threshold_length'] = 1e-1\n",
    "print(eptm.Nf, eptm.Ne, eptm.Nv)\n",
    "draw_spec['face']['color'] = eptm.face_df.area\n",
    "draw_spec['face']['color_range'] = (0, 1)\n",
    "\n",
    "\n",
    "class NormPerimeterElasticity(effectors.PerimeterElasticity):\n",
    "   \n",
    "    @staticmethod\n",
    "    def gradient(eptm):\n",
    "        \n",
    "        grad_srce_, grad_trgt_ = effectors.PerimeterElasticity.gradient(eptm)\n",
    "        ns = to_nd(eptm.upcast_face(eptm.face_df.num_sides).to_numpy(), eptm.dim)\n",
    "        grad_srce, grad_trgt = grad_srce_ * ns/6, grad_trgt_ * ns/6\n",
    "        return grad_srce, grad_trgt\n",
    "    \n",
    "    @staticmethod\n",
    "    def energy(eptm):\n",
    "        \n",
    "        ns = eptm.face_df.num_sides.mean()\n",
    "        return effectors.PerimeterElasticity.energy(eptm) * ns / 6\n",
    "  \n",
    "        \n",
    "\n",
    "effs = [\n",
    "    NormPerimeterElasticity,\n",
    "    effectors.FaceAreaElasticity,\n",
    "    effectors.CellVolumeElasticity\n",
    "]\n",
    "model = model_factory(effs)\n",
    "\n",
    "\n",
    "model_specs =  {\n",
    "    'cell': {\n",
    "        'is_alive': 1,\n",
    "        'vol': 1.0,\n",
    "        'vol_elasticity': 0.1,\n",
    "        'prefered_vol':  6,\n",
    "    },\n",
    "    'face': {\n",
    "        'is_alive': 1,\n",
    "        'area': 1.0,\n",
    "        'area_elasticity': 1.0,\n",
    "        'prefered_area': 1.0,\n",
    "        'perimeter': 1.0,\n",
    "        'perimeter_elasticity': 0.1,\n",
    "        'prefered_perimeter': 3.81\n",
    "    },\n",
    "        'edge': {},\n",
    "        'vert': {},\n",
    "        'settings': {'nrj_norm_factor': 1.0}\n",
    "}\n",
    "\n",
    "eptm.update_specs(model_specs, reset=True)\n",
    "RNRGeometry.update_all(eptm)\n",
    "solver = QSSolver()\n",
    "\n",
    "solver.find_energy_min(eptm, geom, model)\n",
    "\n",
    "ipv.clear()\n",
    "fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I→H transition\n",
    "\n",
    "![IH transition in a bulk and a RNR](data/png/IH_transition.png)\n",
    "\n",
    "This produces a shape similar to the one studied by LM Escudero and collegues [in their article](https://www.nature.com/articles/s41467-018-05376-1.pdf).\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9156c3434b564cf0b9a68bd4a8569cf1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Figure(camera=PerspectiveCamera(fov=45.0, position=(0.0, 0.0, 2.0), projectionMatrix=(1.0, 0.0, 0.0, 0.0, 0.0,…"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_1011 = 300\n",
    "\n",
    "\n",
    "eptm.settings['threshold_length'] = 0.5\n",
    "\n",
    "srce, trgt, face, cell = eptm.edge_df.loc[e_1011,\n",
    "                                          ['srce', 'trgt',\n",
    "                                           'face', 'cell']]\n",
    "\n",
    "IH_transition(eptm, e_1011)\n",
    "RNRGeometry.update_all(eptm)\n",
    "\n",
    "highlight_cells(eptm, cell, reset_visible=True)\n",
    "draw_spec['face']['color'] = eptm.face_df.area\n",
    "\n",
    "_ = solver.find_energy_min(eptm, geom, model)\n",
    "\n",
    "ipv.clear()\n",
    "fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# H → I transition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8b13d8abed784230b569ede821a4ca7f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Figure(camera=PerspectiveCamera(fov=45.0, position=(0.0, 0.0, 2.0), projectionMatrix=(1.0, 0.0, 0.0, 0.0, 0.0,…"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fa = 93\n",
    "\n",
    "\n",
    "HI_transition(eptm, fa)\n",
    "RNRGeometry.update_all(eptm)\n",
    "\n",
    "highlight_cells(eptm, cell, reset_visible=False)\n",
    "ipv.clear()\n",
    "draw_spec['face']['color'] = eptm.face_df.area\n",
    "\n",
    "fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing for I→H / H →I transition triggers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.18569133439530414"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eptm.edge_df.length.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#eptm.settings['threshold_length'] = 1e-3\n",
    "eptm.settings['threshold_length'] = eptm.edge_df.length.min() + 1e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "ih, hi = find_rearangements(eptm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Condition 3\n",
    "\n",
    "This condition is satisfied if `eptm.settings['threshold_length']`\n",
    "is well defined, i.e, small with respect to the unit length or the\n",
    "average edge length.\n",
    "\n",
    "### Condition 4\n",
    "\n",
    "* (i)  Two edges never share two vertices simultaneously.\n",
    "\n",
    "_In our half-edge architecture, this is not as straight forward. But I think the condition can\n",
    "be reformulated as:_ Two edges **from the same face** never share two edges simultanously.\n",
    "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.\n",
    "\n",
    "* (ii) Two polygonal faces never share two or more edges simultaneously.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tyssue.topology.base_topology import condition_4i, condition_4ii"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Int64Index([], dtype='int64', name='face'),\n",
       " array([], shape=(0, 2), dtype=int64))"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "condition_4i(eptm), condition_4ii(eptm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d73f1e6becb4dfa85a31022010efab6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Figure(camera=PerspectiveCamera(fov=45.0, position=(0.0, 0.0, 2.0), projectionMatrix=(1.0, 0.0, 0.0, 0.0, 0.0,…"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "IH_transition(eptm, ih[0])\n",
    "RNRGeometry.update_all(eptm)\n",
    "solver.find_energy_min(eptm, RNRGeometry, model)\n",
    "draw_spec['face']['color'] = eptm.face_df.area\n",
    "\n",
    "ipv.clear()\n",
    "fig, meshes = sheet_view(eptm, mode='3D', **draw_spec)\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert eptm.validate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "widgets": {
   "state": {
    "2071d16c0db447cb9f98bb491bf52bad": {
     "views": [
      {
       "cell_index": 2
      }
     ]
    },
    "4db4a0e7a4f047df8f900b6faa25afe9": {
     "views": [
      {
       "cell_index": 13
      }
     ]
    },
    "57b8348131ba47a2a18f23920a8e2a9b": {
     "views": [
      {
       "cell_index": 6
      }
     ]
    },
    "d9373ebfbb5542c695d896381b0aeefd": {
     "views": [
      {
       "cell_index": 4
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}