Compute curvature on mesh boundary and use in model specifications

examples/example2-3d/example2-3d.ipynb

@@ -128,7 +128,7 @@
     "# Species\n",
     "# =============================================================================================\n",
     "# name, initial concentration, concentration units, diffusion, diffusion units, compartment\n",
-    "A = Species(\"A\", 1.0, vol_unit, 10.0, D_unit, \"Cyto\")\n",
+    "A = Species(\"A\", 1.0, vol_unit, 1.0, D_unit, \"Cyto\")\n",
     "X = Species(\"X\", 1000, surf_unit, 0.1, D_unit, \"PM\")\n",
     "B = Species(\"B\", 0.0, surf_unit, 0.01, D_unit, \"PM\")\n",
     "sc = SpeciesContainer()\n",

examples/example2/example2.ipynb

@@ -518,7 +518,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.6"
+   "version": "3.10.12"
   },
   "vscode": {
    "interpreter": {

examples/example2/example2_withcurv.ipynb

@@ -0,0 +1,377 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "f65f18d7",
+   "metadata": {},
+   "source": [
+    "# Example 2: Simple 2D cell signaling model - with curvature-sensitive reactions\n",
+    "\n",
+    "We model a reaction between the cell interior and cell membrane in a 2D geometry:\n",
+    "- Cyto - 2D cell \"volume\"\n",
+    "- PM - 1D cell boundary (represents plasma membrane)\n",
+    "\n",
+    "We use the model from [Rangamani et al, 2013, Cell](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3874130/), with a slight variation to illustrate curvature-sensitive reactions. A cytosolic species, \"A\", reacts with a species on the PM, \"X\", to form a new species on the PM, \"B\". We consider four different cases of curvature-sensitivity:\n",
+    "1. No curvature sensitivity (reference case)\n",
+    "2. Binding of A to X preferentially occurs in regions of low curvature.\n",
+    "3. Binding of A to X preferentially occurs in regions of high curvature.\n",
+    "4. X is initially localized to regions of high curvature and cannot diffuse.\n",
+    "\n",
+    "In cases 1 and 4, reactions are unaltered from the original model. In cases 2 and 3, reaction definitions include a dependence on curvature. In the original model, A and X binding occurred at rate $k_{on} C_A N_X$, whereas in our altered model, binding occurs at rate $k_{on} \\exp(K \\text{``curv''}) C_A N_X$. $K$ is a constant that controls the curvature sensitivity (negative for low-curvature biased reaction, positive for high-curvature biased reaction) and $\\text{``curv''}$ is the mean curvature computed over the boundary of the mesh.\n",
+    "\n",
+    "The resulting PDE and boundary condition for species A are as follows:\n",
+    "\n",
+    "$$\n",
+    "\\frac{\\partial{C_A}}{\\partial{t}} = D_A \\nabla ^2 C_A \\quad \\text{in} \\; \\Omega_{Cyto}\\\\\n",
+    "\\text{B.C. for A:} \\quad D_A (\\textbf{n} \\cdot \\nabla C_A)  = -k_{on} \\exp(K \\text{``curv''}) C_A N_X + k_{off} N_B \\quad \\text{on} \\; \\Gamma_{PM}\n",
+    "$$\n",
+    "\n",
+    "Similarly, the PDEs for X and B are given by:\n",
+    "$$\n",
+    "\\frac{\\partial{N_X}}{\\partial{t}} = D_X \\nabla ^2 N_X - k_{on} \\exp(K \\text{``curv''}) C_A N_X + k_{off} N_B \\quad \\text{on} \\; \\Gamma_{PM}\\\\\n",
+    "\\frac{\\partial{N_B}}{\\partial{t}} = D_B \\nabla ^2 N_B + k_{on} \\exp(K \\text{``curv''}) C_A N_X - k_{off} N_B \\quad \\text{on} \\; \\Gamma_{PM}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b224bea7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib import pyplot as plt\n",
+    "import matplotlib.image as mpimg\n",
+    "img_A = mpimg.imread('axb-diagram.png')\n",
+    "plt.imshow(img_A)\n",
+    "plt.axis('off')"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "a59f3428",
+   "metadata": {},
+   "source": [
+    "Imports and logger initialization:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cc398816",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import dolfin as d\n",
+    "import sympy as sym\n",
+    "import numpy as np\n",
+    "import pathlib\n",
+    "import logging\n",
+    "import gmsh  # must be imported before pyvista if dolfin is imported first\n",
+    "\n",
+    "from smart import config, common, mesh, model, mesh_tools, visualization\n",
+    "from smart.units import unit\n",
+    "from smart.model_assembly import (\n",
+    "    Compartment,\n",
+    "    Parameter,\n",
+    "    Reaction,\n",
+    "    Species,\n",
+    "    SpeciesContainer,\n",
+    "    ParameterContainer,\n",
+    "    CompartmentContainer,\n",
+    "    ReactionContainer,\n",
+    ")\n",
+    "from matplotlib import pyplot as plt\n",
+    "import matplotlib.image as mpimg\n",
+    "\n",
+    "logger = logging.getLogger(\"smart\")\n",
+    "logger.setLevel(logging.INFO)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "95b9d865",
+   "metadata": {},
+   "source": [
+    "First, we define the various units for use in the model."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4f4023cf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "um = unit.um\n",
+    "molecule = unit.molecule\n",
+    "sec = unit.sec\n",
+    "dimensionless = unit.dimensionless\n",
+    "D_unit = um**2 / sec\n",
+    "surf_unit = molecule / um**2\n",
+    "flux_unit = molecule / (um * sec)\n",
+    "edge_unit = molecule / um"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "46582d26",
+   "metadata": {},
+   "source": [
+    "Next we generate the model by assembling the compartment, species, parameter, and reaction containers (see Example 1 or API documentation for more details)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "09079b17",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# =============================================================================================\n",
+    "# Compartments\n",
+    "# =============================================================================================\n",
+    "# name, topological dimensionality, length scale units, marker value\n",
+    "Cyto = Compartment(\"Cyto\", 2, um, 1)\n",
+    "PM = Compartment(\"PM\", 1, um, 3)\n",
+    "cc = CompartmentContainer()\n",
+    "cc.add([Cyto, PM])\n",
+    "\n",
+    "# =============================================================================================\n",
+    "# Species\n",
+    "# =============================================================================================\n",
+    "# name, initial concentration, concentration units, diffusion, diffusion units, compartment\n",
+    "A = Species(\"A\", 1.0, surf_unit, 10.0, D_unit, \"Cyto\")\n",
+    "X = Species(\"X\", 1.0, edge_unit, 1.0, D_unit, \"PM\")\n",
+    "B = Species(\"B\", 0.0, edge_unit, 1.0, D_unit, \"PM\")\n",
+    "sc = SpeciesContainer()\n",
+    "sc.add([A, X, B])\n",
+    "\n",
+    "# =============================================================================================\n",
+    "# Parameters and Reactions\n",
+    "# =============================================================================================\n",
+    "\n",
+    "# Reaction of A and X to make B (Cyto-PM reaction)\n",
+    "kon = Parameter(\"kon\", 1.0, 1/(surf_unit*sec))\n",
+    "koff = Parameter(\"koff\", 1.0, 1/sec)\n",
+    "r1 = Reaction(\"r1\", [\"A\", \"X\"], [\"B\"],\n",
+    "              param_map={\"on\": \"kon\", \"off\": \"koff\"},\n",
+    "              species_map={\"A\": \"A\", \"X\": \"X\", \"B\": \"B\"},\n",
+    "              eqn_f_str=f\"A*X*on - B*off\")\n",
+    "\n",
+    "pc = ParameterContainer()\n",
+    "pc.add([kon, koff])\n",
+    "rc = ReactionContainer()\n",
+    "rc.add([r1])"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "15c35d39",
+   "metadata": {},
+   "source": [
+    "Now we create a circular mesh (mesh built using gmsh in `smart.mesh_tools`), along with marker functions `mf2` and `mf1`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fe56e162",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create mesh\n",
+    "h_ellipse = 0.1\n",
+    "xrad = 2.0\n",
+    "yrad = 0.5\n",
+    "surf_tag = 1\n",
+    "edge_tag = 3\n",
+    "ellipse_mesh, mf1, mf2 = mesh_tools.create_ellipses(xrad, yrad, hEdge=h_ellipse,\n",
+    "                                                    outer_tag=surf_tag, outer_marker=edge_tag)\n",
+    "# compute curvature at the pm (idenfitied by \"edge_tag\")\n",
+    "mf_curv = mesh_tools.compute_curvature(ellipse_mesh, mf1, mf2, [edge_tag], [surf_tag])\n",
+    "visualization.plot_dolfin_mesh(ellipse_mesh, mf2, view_xy=True)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "fe04ad6b",
+   "metadata": {},
+   "source": [
+    "Write mesh and meshfunctions to file, then create `mesh.ParentMesh` object."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e15255a1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mesh_folder = pathlib.Path(\"ellipse_mesh\")\n",
+    "mesh_folder.mkdir(exist_ok=True)\n",
+    "mesh_file = mesh_folder / \"ellipse_mesh.h5\"\n",
+    "mesh_tools.write_mesh(ellipse_mesh, mf1, mf2, mesh_file)\n",
+    "# save curvatures for reference\n",
+    "curv_file_name = mesh_folder / \"curvatures.xdmf\"\n",
+    "curv_file = d.XDMFFile(str(curv_file_name))\n",
+    "curv_file.write(mf_curv)\n",
+    "\n",
+    "parent_mesh = mesh.ParentMesh(\n",
+    "    mesh_filename=str(mesh_file),\n",
+    "    mesh_filetype=\"hdf5\",\n",
+    "    name=\"parent_mesh\",\n",
+    "    curvature=mf_curv,\n",
+    ")"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "d1c0cab2",
+   "metadata": {},
+   "source": [
+    "Initialize config."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b059df37",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "config_cur = config.Config()\n",
+    "config_cur.solver.update(\n",
+    "    {\n",
+    "        \"final_t\": 5.0,\n",
+    "        \"initial_dt\": 0.05,\n",
+    "        \"time_precision\": 6,\n",
+    "    }\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "adaa03b9",
+   "metadata": {},
+   "source": [
+    "Here, we compare the solution for different cases of curvature sensitivity vs. curvature insensitivity."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "434c7b51",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# folder for saving final mesh images\n",
+    "meshimg_folder = pathlib.Path(\"mesh_images\")\n",
+    "meshimg_folder = meshimg_folder.resolve()\n",
+    "meshimg_folder.mkdir(exist_ok=True)\n",
+    "images = {}\n",
+    "\n",
+    "curv_const = [0.0, -0.5, 0.5, np.nan]\n",
+    "label_str = [\"No curvature sensitivity\", \"Low curvature reaction\", \n",
+    "             \"High curvature reaction\", \"Curvature-dependent distribution of X\"]\n",
+    "for i in range(len(curv_const)):\n",
+    "    if np.isnan(curv_const[i]):\n",
+    "        r1.eqn_f_str = f\"A*X*on - B*off\"\n",
+    "        X.initial_condition = \"exp(0.5*curv)\"\n",
+    "        X.initial_condition_expression = \"exp(0.5*curv)\"\n",
+    "        X.D = 0.0\n",
+    "    else:\n",
+    "        r1.eqn_f_str = f\"A*X*exp({curv_const[i]}*curv)*on - B*off\"\n",
+    "    model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh)\n",
+    "    model_cur.initialize()\n",
+    "    results = dict()\n",
+    "    result_folder = pathlib.Path(f\"resultsEllipse{i}\")\n",
+    "    result_folder.mkdir(exist_ok=True)\n",
+    "    for species_name, species in model_cur.sc.items:\n",
+    "        results[species_name] = d.XDMFFile(\n",
+    "            model_cur.mpi_comm_world, str(result_folder / f\"{species_name}.xdmf\")\n",
+    "        )\n",
+    "        results[species_name].parameters[\"flush_output\"] = True\n",
+    "        results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+    "    avg_A = [A.initial_condition]\n",
+    "    dx = d.Measure(\"dx\", domain=model_cur.cc['Cyto'].dolfin_mesh)\n",
+    "    volume = d.assemble(1.0*dx)\n",
+    "    while True:\n",
+    "        # Solve the system\n",
+    "        model_cur.monolithic_solve()\n",
+    "        # Save results for post processing\n",
+    "        for species_name, species in model_cur.sc.items:\n",
+    "            results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+    "        int_val = d.assemble(model_cur.sc['A'].u['u']*dx)\n",
+    "        avg_A.append(int_val / volume)\n",
+    "        # End if we've passed the final time\n",
+    "        if model_cur.t >= model_cur.final_t:\n",
+    "            break\n",
+    "    plt.plot(model_cur.tvec, avg_A, label=label_str[i])\n",
+    "    meshimg_file = meshimg_folder / f\"ellipse_mesh{i}.png\"\n",
+    "    images[i] = visualization.plot(model_cur.sc[\"A\"].u[\"u\"], view_xy=True, filename=meshimg_file)\n",
+    "plt.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "28c08fa6",
+   "metadata": {},
+   "source": [
+    "Show final condition for all cases."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c5b31511",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(2, 2, figsize=(20, 15))\n",
+    "for axi, (idx, image) in zip(ax.flatten(), images.items()):\n",
+    "    axi.imshow(image)\n",
+    "    axi.axis('off')\n",
+    "    axi.set_title(label_str[idx])"
+   ]
+  }
+ ],
+ "metadata": {
+  "jupytext": {
+   "cell_metadata_filter": "-all",
+   "main_language": "python",
+   "notebook_metadata_filter": "-all"
+  },
+  "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.12"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}