From 1768a6371f6ce729f3046cb4ff8f15be4ea5a98b Mon Sep 17 00:00:00 2001 From: emmetfrancis <99422170+emmetfrancis@users.noreply.github.com> Date: Mon, 23 Jun 2025 11:09:43 -0700 Subject: [PATCH 1/8] Allow use of additional keys to store information about subdomains in mesh, include 2d mesh creation in xy plane (removing z coord in gmsh) --- smart/mesh.py | 37 +++++- smart/mesh_tools.py | 287 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 322 insertions(+), 2 deletions(-) diff --git a/smart/mesh.py b/smart/mesh.py index 96774ce..ba93abd 100644 --- a/smart/mesh.py +++ b/smart/mesh.py @@ -152,6 +152,10 @@ class ParentMesh(_Mesh): as an xml or xdmf file containing a vertex mesh function of type double extra_keys (list): list of names of extra keys to load from hdf5 mesh file, specifying other subdomains besides main compartments (optional) + key_for_cells (string): a string specifying the extra key to use as a + mesh function to define cell markers in mesh (optional) + key_for_facets (string): a string specifying the extra key to use as a + mesh function to define facet markers in mesh (optional) """ mesh_filename: str @@ -161,6 +165,8 @@ class ParentMesh(_Mesh): use_partition: bool curvature: d.MeshFunction extra_keys: list = [] + key_for_cells: str = "" + key_for_facets: str = "" def __init__( self, @@ -171,6 +177,8 @@ def __init__( mpi_comm=d.MPI.comm_world, curvature=None, extra_keys=[], + key_for_cells="", + key_for_facets="", ): super().__init__(name) self.use_partition = use_partition @@ -202,6 +210,30 @@ def __init__( # Otherwise just take what we got self.curvature = curvature + # set cells and/or facets according to extra keys if applicable + if key_for_cells != "": + try: + idx = self.extra_keys.index(key_for_cells) + except ValueError: + raise ValueError(f"'{key_for_cells}' does not match an extra key") + assert self.subdomains[idx].dim() == self.dimensionality, ( + f"Mesh function associated with '{key_for_cells}' " + "does not match mesh cell dimension" + ) + self.mf["cells"] = self.subdomains[idx] + logger.info(f"Cell mesh function loaded from key '{key_for_cells}'") + if key_for_facets != "": + try: + idx = self.extra_keys.index(key_for_facets) + except ValueError: + raise ValueError(f"'{key_for_cells}' does not match an extra key") + assert self.subdomains[idx].dim() == self.dimensionality - 1, ( + f"Mesh function associated with '{key_for_facets}' " + "does not match mesh facet dimension" + ) + self.mf["facets"] = self.subdomains[idx] + logger.info(f"Facet mesh function loaded from key '{key_for_facets}'") + def get_mesh_from_id(self, id): "Find the mesh that has the matching id." # find the mesh in that has the matching id @@ -290,8 +322,9 @@ def read_parent_mesh_functions_from_file(self): assert len(self.child_meshes) > 0 # Init mesh functions - self.mf["cells"] = self._read_parent_mesh_function_from_file(volume_dim) - if self.has_surface: + if "cells" not in self.mf.keys(): + self.mf["cells"] = self._read_parent_mesh_function_from_file(volume_dim) + if self.has_surface and "facets" not in self.mf.keys(): self.mf["facets"] = self._read_parent_mesh_function_from_file(surface_dim) # If any cell markers are given as a list we also create mesh diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py index e1a700a..fb04478 100644 --- a/smart/mesh_tools.py +++ b/smart/mesh_tools.py @@ -1300,6 +1300,291 @@ def meshSizeCallback(dim, tag, x, y, z, lc): return (dmesh, mf2, mf3) +def create_2Dcell_xy( + outerExpr: str = "", + innerExpr: str = "", + hEdge: float = 0, + hInnerEdge: float = 0, + interface_marker: int = 12, + outer_marker: int = 10, + inner_tag: int = 2, + outer_tag: int = 1, + comm: MPI.Comm = d.MPI.comm_world, + verbose: bool = False, + half_cell: bool = True, + return_curvature: bool = False, +) -> Tuple[d.Mesh, d.MeshFunction, d.MeshFunction]: + """ + Creates a 2D mesh of a cell profile, with the bounding curve defined in + terms of r and z (e.g. unit circle would be "r**2 + (z-1)**2 - 1) + It is assumed that substrate is present at z = 0, so if the curve extends + below z = 0 , there is a sharp cutoff. + If half_cell = True, only have of the contour is constructed, with a + left zero-flux boundary at r = 0. + Can include one compartment inside another compartment. + Recommended for use with the axisymmetric feature of SMART. + + Args: + outerExpr: String implicitly defining an r-z curve for the outer surface + innerExpr: String implicitly defining an r-z curve for the inner surface + hEdge: maximum mesh size at the outer edge + hInnerEdge: maximum mesh size at the edge + of the inner compartment + interface_marker: The value to mark facets on the interface with + outer_marker: The value to mark facets on edge of the outer ellipse with + inner_tag: The value to mark the inner ellipse surface with + outer_tag: The value to mark the outer ellipse surface with + comm: MPI communicator to create the mesh with + verbose: If true print gmsh output, else skip + half_cell: If true, consider r=0 the symmetry axis for an axisymm shape + Returns: + Tuple (mesh, facet_marker, cell_marker) + """ + import gmsh + + if outerExpr == "": + ValueError("Outer surface is not defined") + + if return_curvature: + # create full mesh for curvature analysis and then map onto half mesh + # if half_cell_with_curvature is True + half_cell_with_curvature = half_cell + half_cell = False + + rValsOuter, zValsOuter = implicit_curve(outerExpr) + + if not innerExpr == "": + rValsInner, zValsInner = implicit_curve(innerExpr) + zMid = np.mean(zValsInner) + ROuterVec = np.sqrt(rValsOuter**2 + (zValsOuter - zMid) ** 2) + RInnerVec = np.sqrt(rValsInner**2 + (zValsInner - zMid) ** 2) + maxOuterDim = max(ROuterVec) + maxInnerDim = max(RInnerVec) + else: + zMid = np.mean(zValsOuter) + ROuterVec = np.sqrt(rValsOuter**2 + (zValsOuter - zMid) ** 2) + maxOuterDim = max(ROuterVec) + if np.isclose(hEdge, 0): + hEdge = 0.1 * maxOuterDim + if np.isclose(hInnerEdge, 0): + hInnerEdge = 0.2 * maxOuterDim if innerExpr == "" else 0.2 * maxInnerDim + # Create the 2D mesh using gmsh + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", int(verbose)) + gmsh.model.add("2DCell") + # first add outer body and revolve + outer_tag_list = [] + for i in range(len(rValsOuter)): + cur_tag = gmsh.model.occ.add_point(rValsOuter[i], zValsOuter[i], 0.0) + outer_tag_list.append(cur_tag) + outer_spline_tag = gmsh.model.occ.add_spline(outer_tag_list) + if not half_cell: + outer_tag_list2 = [] + for i in range(len(rValsOuter)): + cur_tag = gmsh.model.occ.add_point(-rValsOuter[i], zValsOuter[i], 0.0) + outer_tag_list2.append(cur_tag) + outer_spline_tag2 = gmsh.model.occ.add_spline(outer_tag_list2) + if np.isclose(zValsOuter[-1], 0): # then include substrate at z=0 + if half_cell: + origin_tag = gmsh.model.occ.add_point(0, 0, 0) + symm_axis_tag = gmsh.model.occ.add_line(origin_tag, outer_tag_list[0]) + bottom_tag = gmsh.model.occ.add_line(origin_tag, outer_tag_list[-1]) + outer_loop_tag = gmsh.model.occ.add_curve_loop( + [outer_spline_tag, bottom_tag, symm_axis_tag] + ) + else: + bottom_tag = gmsh.model.occ.add_line(outer_tag_list[-1], outer_tag_list2[-1]) + outer_loop_tag = gmsh.model.occ.add_curve_loop( + [outer_spline_tag, outer_spline_tag2, bottom_tag] + ) + else: + if half_cell: + symm_axis_tag = gmsh.model.occ.add_line(outer_tag_list[0], outer_tag_list[-1]) + outer_loop_tag = gmsh.model.occ.add_curve_loop([outer_spline_tag, symm_axis_tag]) + else: + outer_loop_tag = gmsh.model.occ.add_curve_loop([outer_spline_tag, outer_spline_tag2]) + cell_plane_tag = gmsh.model.occ.add_plane_surface([outer_loop_tag]) + + if innerExpr == "": + # No inner shape in this case + gmsh.model.occ.synchronize() + gmsh.model.add_physical_group(2, [cell_plane_tag], tag=outer_tag) + facets = gmsh.model.getBoundary([(2, cell_plane_tag)]) + facet_tag_list = [] + for i in range(len(facets)): + facet_tag_list.append(facets[i][1]) + if half_cell: # if half, set symmetry axis to 0 (no flux) + xmin, ymin, zmin = (-hInnerEdge / 10, -1, -hInnerEdge / 10) + xmax, ymax, zmax = (hInnerEdge / 10, max(zValsOuter) + 1, hInnerEdge / 10) + all_symm_bound = gmsh.model.occ.get_entities_in_bounding_box( + xmin, ymin, zmin, xmax, ymax, zmax, dim=1 + ) + symm_bound_markers = [] + for i in range(len(all_symm_bound)): + symm_bound_markers.append(all_symm_bound[i][1]) + gmsh.model.add_physical_group(1, symm_bound_markers, tag=0) + gmsh.model.add_physical_group(1, facet_tag_list, tag=outer_marker) + else: + # Add inner shape + inner_tag_list = [] + for i in range(len(rValsInner)): + cur_tag = gmsh.model.occ.add_point(rValsInner[i], zValsInner[i], 0.0) + inner_tag_list.append(cur_tag) + inner_spline_tag = gmsh.model.occ.add_spline(inner_tag_list) + if half_cell: + symm_inner_tag = gmsh.model.occ.add_line(inner_tag_list[0], inner_tag_list[-1]) + inner_loop_tag = gmsh.model.occ.add_curve_loop([inner_spline_tag, symm_inner_tag]) + else: + inner_tag_list2 = [] + for i in range(len(rValsInner)): + cur_tag = gmsh.model.occ.add_point(-rValsInner[i], zValsInner[i], 0.0) + inner_tag_list2.append(cur_tag) + inner_spline_tag2 = gmsh.model.occ.add_spline(inner_tag_list2) + inner_loop_tag = gmsh.model.occ.add_curve_loop([inner_spline_tag, inner_spline_tag2]) + inner_plane_tag = gmsh.model.occ.add_plane_surface([inner_loop_tag]) + cell_plane_list = [cell_plane_tag] + inner_plane_list = [inner_plane_tag] + + outer_volume = [] + inner_volume = [] + all_volumes = [] + inner_marker_list = [] + outer_marker_list = [] + for i in range(len(cell_plane_list)): + cell_plane_tag = cell_plane_list[i] + inner_plane_tag = inner_plane_list[i] + # Create interface between 2 objects + two_shapes, (outer_shape_map, inner_shape_map) = gmsh.model.occ.fragment( + [(2, cell_plane_tag)], [(2, inner_plane_tag)] + ) + gmsh.model.occ.synchronize() + + # Get the outer boundary + outer_shell = gmsh.model.getBoundary(two_shapes, oriented=False) + for i in range(len(outer_shell)): + outer_marker_list.append(outer_shell[i][1]) + # Get the inner boundary + inner_shell = gmsh.model.getBoundary(inner_shape_map, oriented=False) + for i in range(len(inner_shell)): + inner_marker_list.append(inner_shell[i][1]) + for tag in outer_shape_map: + all_volumes.append(tag[1]) + for tag in inner_shape_map: + inner_volume.append(tag[1]) + + for vol in all_volumes: + if vol not in inner_volume: + outer_volume.append(vol) + + # Add physical markers for facets + if half_cell: # if half, set symmetry axis to 0 (no flux) + xmin, ymin, zmin = (-hInnerEdge / 10, -1, -hInnerEdge / 10) + xmax, ymax, zmax = (hInnerEdge / 10, max(zValsOuter) + 1, hInnerEdge / 10) + all_symm_bound = gmsh.model.occ.get_entities_in_bounding_box( + xmin, ymin, zmin, xmax, ymax, zmax, dim=1 + ) + symm_bound_markers = [] + for i in range(len(all_symm_bound)): + symm_bound_markers.append(all_symm_bound[i][1]) + # note that this first call sets the symmetry axis to tag 0 and + # this is not overwritten by the next calls to add_physical_group + gmsh.model.add_physical_group(1, symm_bound_markers, tag=0) + gmsh.model.add_physical_group(1, outer_marker_list, tag=outer_marker) + gmsh.model.add_physical_group(1, inner_marker_list, tag=interface_marker) + + # Physical markers for "volumes" + gmsh.model.add_physical_group(2, outer_volume, tag=outer_tag) + gmsh.model.add_physical_group(2, inner_volume, tag=inner_tag) + + def meshSizeCallback(dim, tag, x, y, z, lc): + # mesh length is hEdge at the PM and hInnerEdge at the inner membrane + # between these, the value is interpolated based on the relative distance + # between the two membranes. + # Inside the inner shape, the value is interpolated between hInnerEdge + # and lc3, where lc3 = max(hInnerEdge, 0.2*maxInnerDim) + # if innerRad=0, then the mesh length is interpolated between + # hEdge at the PM and 0.2*maxOuterDim in the center + rCur = np.sqrt(x**2 + z**2) + RCur = np.sqrt(rCur**2 + (y - zMid) ** 2) + outer_dist = np.sqrt((rCur - rValsOuter) ** 2 + (z - zValsOuter) ** 2) + np.append(outer_dist, z) # include the distance from the substrate + dist_to_outer = min(outer_dist) + if innerExpr == "": + lc3 = 0.2 * maxOuterDim + dist_to_inner = RCur + in_outer = True + else: + inner_dist = np.sqrt((rCur - rValsInner) ** 2 + (y - zValsInner) ** 2) + dist_to_inner = min(inner_dist) + inner_idx = np.argmin(inner_dist) + inner_rad = RInnerVec[inner_idx] + R_rel_inner = RCur / inner_rad + lc3 = max(hInnerEdge, 0.2 * maxInnerDim) + in_outer = R_rel_inner > 1 + lc1 = hEdge + lc2 = hInnerEdge + if in_outer: + lcTest = lc1 + (lc2 - lc1) * (dist_to_outer) / (dist_to_inner + dist_to_outer) + else: + lcTest = lc2 + (lc3 - lc2) * (1 - R_rel_inner) + return lcTest + + gmsh.model.mesh.setSizeCallback(meshSizeCallback) + # set off the other options for mesh size determination + gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + # this changes the algorithm from Frontal-Delaunay to Delaunay, + # which may provide better results when there are larger gradients in mesh size + gmsh.option.setNumber("Mesh.Algorithm", 5) + + gmsh.model.mesh.generate(2) + rank = MPI.COMM_WORLD.rank + tmp_folder = pathlib.Path(f"tmp_2DCell_{rank}") + tmp_folder.mkdir(exist_ok=True) + gmsh_file = tmp_folder / "2DCell.msh" + gmsh.write(str(gmsh_file)) + gmsh.finalize() + + # return dolfin mesh of max dimension (parent mesh) and marker functions mf2 and mf3 + dmesh, mf2, mf3 = gmsh_to_dolfin(str(gmsh_file), tmp_folder, 2, comm) + # remove tmp mesh and tmp folder + gmsh_file.unlink(missing_ok=False) + tmp_folder.rmdir() + # return dolfin mesh, mf2 (2d tags) and mf3 (3d tags) + if return_curvature: + if innerExpr == "": + facet_list = [outer_marker] + cell_list = [outer_tag] + else: + facet_list = [outer_marker, interface_marker] + cell_list = [outer_tag, inner_tag] + if half_cell_with_curvature: # will likely not work in parallel... + dmesh_half, mf2_half, mf3_half = create_2Dcell( + outerExpr, + innerExpr, + hEdge, + hInnerEdge, + interface_marker, + outer_marker, + inner_tag, + outer_tag, + comm, + verbose, + half_cell=True, + return_curvature=False, + ) + kappa_mf = compute_curvature( + dmesh, mf2, mf3, facet_list, cell_list, half_mesh_data=(dmesh_half, mf2_half) + ) + (dmesh, mf2, mf3) = (dmesh_half, mf2_half, mf3_half) + else: + kappa_mf = compute_curvature(dmesh, mf2, mf3, facet_list, cell_list) + return (dmesh, mf2, mf3, kappa_mf) + else: + return (dmesh, mf2, mf3) + + def gmsh_to_dolfin( gmsh_file_name: str, tmp_folder: pathlib.Path = pathlib.Path("tmp_folder"), @@ -1340,6 +1625,8 @@ def gmsh_to_dolfin( # convert cell mesh cells = mesh_in.get_cells_type(cell_type) cell_data = mesh_in.get_cell_data("gmsh:physical", cell_type) # extract values of tags + if dimension == 2 and np.all(mesh_in.points[:, 2] == 0): + mesh_in.points = mesh_in.points[:, :2] # then prune z values out_mesh_cell = meshio.Mesh( points=mesh_in.points, cells={cell_type: cells}, From 10aef259a3e34919fa37255fa722930fdbeeb634 Mon Sep 17 00:00:00 2001 From: emmetfrancis Date: Mon, 23 Jun 2025 22:30:24 -0700 Subject: [PATCH 2/8] New feature to allow parameters with assigned values over a region of the mesh, introduce new version of example 5 to illustrate SOCE --- examples/example5/example5_withSOCE.ipynb | 837 ++++++++++++++++++++++ examples/example5/model_cur.pkl | Bin 0 -> 3487 bytes smart/model.py | 10 + smart/model_assembly.py | 34 +- 4 files changed, 880 insertions(+), 1 deletion(-) create mode 100644 examples/example5/example5_withSOCE.ipynb create mode 100644 examples/example5/model_cur.pkl diff --git a/examples/example5/example5_withSOCE.ipynb b/examples/example5/example5_withSOCE.ipynb new file mode 100644 index 0000000..d1405bd --- /dev/null +++ b/examples/example5/example5_withSOCE.ipynb @@ -0,0 +1,837 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "f65f18d7", + "metadata": {}, + "source": [ + "# Example 5: Generic cell signaling system in 3D, illustrating a SOCE-like effect\n", + "\n", + "This example is slightly altered from the normal version of example 5 to illustrate the role of store-operated calcium entry (SOCE) accompanying calcium release from an interior compartment.\n", + "\n", + "Geometry is divided into 4 domains; two volumes, and two surfaces:\n", + "- cytosol (Cyto): $\\Omega_{Cyto}$\n", + "- endoplasmic reticulum volume (ER): $\\Omega_{ER}$\n", + "- plasma membrane (PM): $\\Gamma_{PM}$\n", + "- ER membrane (ERm): $\\Gamma_{ERm}$\n", + "\n", + "For simplicity, here we consider a \"cube-within-a-cube\" geometry, in which the smaller\n", + "inner cube represents a section of ER and one face of the outer cube ($x=0$) represents the PM. The other\n", + "faces of the outer cube are treated as no flux boundaries. The space outside\n", + "the inner cube but inside the outer cube is classified as cytosol.\n", + "\n", + "There are three function-spaces on these three domains:\n", + "\n", + "$$\n", + "u^{Cyto} = [A, B] \\quad \\text{on} \\quad \\Omega^{Cyto}\\\\\n", + "u^{ER} = [AER] \\quad \\text{on} \\quad \\Omega^{ER}\\\\\n", + "v^{ERm} = [R, Ro] \\quad \\text{on} \\quad \\Gamma^{ERm}\n", + "$$\n", + "\n", + "In words, this says that species A and B reside in the cytosolic volume, \n", + "species AER corresponds to an amount of species A that lives in the ER volume,\n", + "and species R (closed receptor/channel) and Ro (open receptor/channel) reside on the ER membrane.\n", + "\n", + "In this model, species B reacts with a receptor/channel, R, on the ER membrane, causing it to open (change state from R->Ro), \n", + "allowing species A to flow out of the ER and into the cytosol. \n", + "Note that this is roughly similar to an IP3 pulse at the PM, leading to Ca2+ release from the ER,\n", + "where, by analogy, species B is similar to IP3 and species A is similar to Ca2+. A more comprehensive\n", + "model of Ca2+ dynamics in particular is implemented in Example 6.\n", + "\n", + "Here, we also introduce a flux analogous to SOCE, wherein the flux of A through the PM depends on depletion of A from the ER.\n", + "For simplicity, we scale this flux with the distance between the PM and ER ($d_{PM-ER}$) and the concentration of A at closest point of ER to the PM ($A_{ER,close}$):\n", + "\n", + "$$\n", + "J_{SOCE} = \\frac{d_0}{d_{PM-ER} + d_0} \\frac{1}{1 + \\left(\\frac{A_{ER}}{A_{ref}}\\right)^4},\n", + "$$\n", + "\n", + "where $d_0$ is a characteristic length scale controlling the likelihood of STIM-Orai binding and $A_{ref}$ governs the ER calcium depletion required to trigger SOCE." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5edff5d", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "import matplotlib.image as mpimg\n", + "img_A = mpimg.imread('example5-diagram.png')\n", + "plt.imshow(img_A)\n", + "plt.axis('off')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "f543d75e", + "metadata": {}, + "source": [ + "As specified in our [mathematical documentation](https://rangamanilabucsd.github.io/smart/docs/math.html), assuming diffusive transport, the PDE and boundary condition for each of these volumetric species takes the form:\n", + "\n", + "$$\n", + "\\partial_t u_i^m - \\nabla \\cdot ( D_i^m \\nabla (u_i^m) ) - f_i^m(u_i^m) = 0 \\qquad \\text{ in } \\Omega^m\\\\\n", + "D_i \\nabla u_i^m \\cdot n^m - R_i^{q} (u^m, u^n, v^q) = 0 \\qquad \\text{ on } \\Gamma^{q}\n", + "$$\n", + "\n", + "and the surface species take the form:\n", + "\n", + "$$\n", + "\\partial_t v_i^q - \\nabla_S \\cdot (D_i^q \\nabla_S v ) - g_i^q ( u^m, u^n, v^q ) = 0 \\qquad \\text{ on } \\Gamma^{q}\\\\\n", + "D_i \\nabla v_i^q \\cdot n^q = 0 \\qquad \\text{ on } \\partial\\Gamma^{q}\n", + "$$\n", + "\n", + "Our reaction terms and boundary conditions are chosen according to the system described above. For the purposes of this simplified example we use linear mass action in all reaction terms except SOCE. Explicitly writing out this system of PDEs, we have:\n", + "\n", + "\\begin{align}\n", + " \\partial_t u_B^{Cyto} - D_B^{Cyto} \\nabla^2 u_B^{Cyto} + k_{2f} u_B^{Cyto} &= 0 \\qquad \\text{ in } \\Omega^{Cyto}\\\\\n", + " D_B^{Cyto} \\nabla u_B^{Cyto} \\cdot n^{Cyto} + j_1[t] &= 0 \\qquad \\text{ on } \\Gamma^{PM} \\nonumber\\\\\n", + " D_B^{Cyto} \\nabla u_B^{Cyto} \\cdot n^{Cyto} + k_{3f} v_R^{ERm} u_B^{Cyto} - k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n", + " \\nonumber \\\\\n", + " \\partial_t u_A^{Cyto} - D_A^{Cyto} \\nabla^2 u_A^{Cyto} &= 0 \\qquad \\text{ in } \\Omega^{Cyto}\\\\\n", + " D_A^{Cyto} \\nabla u_A^{Cyto} \\cdot n^{Cyto} - J_{SOCE} &= 0 \\qquad \\text{ on } \\Gamma^{PM} \\nonumber\\\\\n", + " D_A^{Cyto} \\nabla u_A^{Cyto} \\cdot n^{Cyto} - k_{4,Vmax} v_{Ro}^{ERm} (u_{AER}^{ER} - u_A^{Cyto}) &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n", + " \\nonumber \\\\\n", + " \\partial_t u_{AER}^{ER} - D_{AER}^{ER} \\nabla^2 u_{AER}^{ER} &= 0 \\qquad \\text{ in } \\Omega^{ER}\\\\\n", + " D_{AER}^{ER} \\nabla u_{AER}^{ER} \\cdot n^{ER} + k_{4,Vmax} v_{Ro}^{ERm} (u_{AER}^{ER} - u_A^{Cyto}) &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n", + " \\nonumber \\\\\n", + " \\partial_t v_{R}^{ERm} - D_{R}^{ERm} \\nabla^2 v_{R}^{ERm} - \n", + " k_{3f} v_R^{ERm} u_B^{Cyto} + k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm}\\\\\n", + " \\nonumber \\\\\n", + " \\partial_t v_{Ro}^{ERm} - D_{Ro}^{ERm} \\nabla^2 v_{Ro}^{ERm} +\n", + " k_{3f} v_R^{ERm} u_B^{Cyto} - k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm}\\\\\n", + "\\end{align}" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "28f77cbf", + "metadata": {}, + "source": [ + "## Code imports and initialization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc398816", + "metadata": {}, + "outputs": [], + "source": [ + "import logging\n", + "\n", + "import dolfin as d\n", + "import sympy as sym\n", + "import numpy as np\n", + "import pathlib\n", + "import gmsh # must be imported before pyvista if dolfin is imported first\n", + "\n", + "from smart import config, mesh, model, mesh_tools, visualization\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 smart.units import unit" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c8650536", + "metadata": {}, + "source": [ + "We will set the logging level to `INFO`. This will display some output during the simulation. If you want to get even more output you could set the logging level to `DEBUG`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ed0899a", + "metadata": {}, + "outputs": [], + "source": [ + "logger = logging.getLogger(\"smart\")\n", + "logger.setLevel(logging.INFO)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "defc1095", + "metadata": {}, + "source": [ + "Futhermore, you could also save the logs to a file by attaching a file handler to the logger as follows.\n", + "\n", + "```\n", + "file_handler = logging.FileHandler(\"filename.log\")\n", + "file_handler.setFormatter(logging.Formatter(smart.config.base_format))\n", + "logger.addHandler(file_handler)\n", + "```" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "95b9d865", + "metadata": {}, + "source": [ + "First, we define the various units for the inputs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f4023cf", + "metadata": {}, + "outputs": [], + "source": [ + "# Aliases - base units\n", + "uM = unit.uM\n", + "um = unit.um\n", + "molecule = unit.molecule\n", + "sec = unit.sec\n", + "# Aliases - units used in model\n", + "D_unit = um**2 / sec\n", + "flux_unit = molecule / (um**2 * sec)\n", + "vol_unit = uM\n", + "surf_unit = molecule / um**2" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "46582d26", + "metadata": {}, + "source": [ + "## Generate model\n", + "Next we generate the model described in the equations above.\n", + "\n", + "### Compartments\n", + "As described above, the three compartments are the cytosol (\"Cyto\"), the plasma membrane (\"PM\"), the ER membrane (\"ERm\"), and the ER interior volume (\"ER\").\n", + "\n", + "Note that, as shown, we can also specify nonadjacency for compartments; this is not strictly necessary, but will generally speed up the simulations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2f34b3c", + "metadata": {}, + "outputs": [], + "source": [ + "Cyto = Compartment(\"Cyto\", 3, um, 1)\n", + "PM = Compartment(\"PM\", 2, um, 10)\n", + "ER = Compartment(\"ER\", 3, um, 2)\n", + "ERm = Compartment(\"ERm\", 2, um, 12)\n", + "PM.specify_nonadjacency(['ERm', 'ER'])\n", + "ERm.specify_nonadjacency(['PM'])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d93ce862", + "metadata": {}, + "source": [ + "Initialize a compartment container and add the 4 compartments to it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "701577e8", + "metadata": {}, + "outputs": [], + "source": [ + "cc = CompartmentContainer()\n", + "cc.add([ERm, ER, PM, Cyto])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "0a6acf0b", + "metadata": {}, + "source": [ + "### Species\n", + "In this case, we have 5 species across 3 different compartments. We define each in turn:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d991278", + "metadata": {}, + "outputs": [], + "source": [ + "A = Species(\"A\", 0.01, vol_unit, 1.0, D_unit, \"Cyto\")\n", + "B = Species(\"B\", 0.0, vol_unit, 1.0, D_unit, \"Cyto\")\n", + "AER = Species(\"AER\", 200.0, vol_unit, 5.0, D_unit, \"ER\")\n", + "# Uniform initial condition of R\n", + "R1 = Species(\"R1\", 1.0, surf_unit, 0.02, D_unit, \"ERm\")\n", + "R1o = Species(\"R1o\", 0.0, surf_unit, 0.02, D_unit, \"ERm\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "b60826cc", + "metadata": {}, + "source": [ + "Create species container and add the 5 species objects to it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e86bebaf", + "metadata": {}, + "outputs": [], + "source": [ + "sc = SpeciesContainer()\n", + "sc.add([R1o, R1, AER, B, A])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "af900a73", + "metadata": {}, + "source": [ + "### Parameters and Reactions\n", + "\n", + "Parameters and reactions are generally defined together, although the order does not strictly matter. We define them in turn as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2120ac6", + "metadata": {}, + "outputs": [], + "source": [ + "# Degradation of B in the cytosol\n", + "k2f = Parameter(\"k2f\", 10, 1 / sec)\n", + "r2 = Reaction(\n", + " \"r2\", [\"B\"], [], param_map={\"on\": \"k2f\"}, reaction_type=\"mass_action_forward\"\n", + ")\n", + "\n", + "# Activating receptors on ERm with B\n", + "k3f = Parameter(\"k3f\", 100, 1 / (uM * sec))\n", + "k3r = Parameter(\"k3r\", 100, 1 / sec)\n", + "r3 = Reaction(\"r3\", [\"B\", \"R1\"], [\"R1o\"], {\"on\": \"k3f\", \"off\": \"k3r\"})\n", + "# Release of A from ERm to cytosol\n", + "k4Vmax = Parameter(\"k4Vmax\", 2000, 1 / (uM * sec))\n", + "r4 = Reaction(\n", + " \"r4\",\n", + " [\"AER\"],\n", + " [\"A\"],\n", + " param_map={\"Vmax\": \"k4Vmax\"},\n", + " species_map={\"R1o\": \"R1o\", \"uER\": \"AER\", \"u\": \"A\"},\n", + " eqn_f_str=\"Vmax*R1o*(uER-u)\",\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "7fdd58b8", + "metadata": {}, + "source": [ + "We define an additional reaction as the time-dependent production of species B at the plasma membrane. In this case, we define a pulse-type function as the derivative of an arctan function. Note that this is useful because we can provide an expression to use for pre-integration.\n", + "\n", + "$$\n", + "j_{int}[t] = V_{max} \\arctan\\left({m (t - t_0)}\\right)\\\\\n", + "j_1[t] = \\frac{dj_{int}[t]}{dt} = \\frac{m V_{max}}{1 + m^2 (t-t_0)^2}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51abf270", + "metadata": {}, + "outputs": [], + "source": [ + "Vmax, t0, m = 500, 0.1, 200\n", + "t = sym.symbols(\"t\")\n", + "pulseI = Vmax * sym.atan(m * (t - t0))\n", + "pulse = sym.diff(pulseI, t)\n", + "j1pulse = Parameter.from_expression(\n", + " \"j1pulse\", pulse, flux_unit, use_preintegration=True, preint_sym_expr=pulseI\n", + ")\n", + "r1 = Reaction(\n", + " \"r1\",\n", + " [],\n", + " [\"B\"],\n", + " param_map={\"J\": \"j1pulse\"},\n", + " eqn_f_str=\"J\",\n", + " explicit_restriction_to_domain=\"PM\",\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "b5e218b0", + "metadata": {}, + "source": [ + "We can plot the time-dependent input by converting the sympy expression to a numpy function using lambdify." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a0829a6b", + "metadata": {}, + "outputs": [], + "source": [ + "from sympy.utilities.lambdify import lambdify\n", + "pulse_func = lambdify(t, pulse, 'numpy') # returns a numpy-ready function\n", + "tArray = np.linspace(0, 1, 100)\n", + "pulse_vals = pulse_func(tArray)\n", + "plt.plot(tArray, pulse_vals)" + ] + }, + { + "cell_type": "markdown", + "id": "55c27105", + "metadata": {}, + "source": [ + "Finally, we define the store-operated calcium (or here, A) influx. \n", + "This is not immediately straightforward, as a flux at the PM depends on concentration of A in the ER.\n", + "One option would be to explicitly introduce contacts between the ER and PM, but this would imply direct entry of calcium into the ER through Orai1.\n", + "Instead, we assume the likelihood of STIM-Orai1 binding scales with PM-ER distance (if they are very close, membrane fluctuations are likely to allow contact between the ERM and PM occasionally, for instance).\n", + "To adopt this strategy, we must construct a map between each point on the PM and the closest point on the ER, which we do after initializing the model below.\n", + "For now, we define the flux as a parameter varying over the PM mesh by calling `Parameter.mesh_quantity()`\n", + "" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a56bbccb", + "metadata": {}, + "outputs": [], + "source": [ + "PSOCE = Parameter.mesh_quantity(\"PSOCE\", 0, unit.dimensionless, compartment=\"PM\")\n", + "J0_SOCE = Parameter(\"J0_SOCE\", 1e4, flux_unit)\n", + "rSOCE = Reaction(\n", + " \"rSOCE\",\n", + " [],\n", + " [\"A\"],\n", + " param_map={\"P\": \"PSOCE\", \"J0\": \"J0_SOCE\"},\n", + " eqn_f_str=\"P*J0\",\n", + " explicit_restriction_to_domain=\"PM\",\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "3d883c6e", + "metadata": {}, + "source": [ + "Create containers for parameters and reactions and add all the parameters and reaction objects to them." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ae2c2c1", + "metadata": {}, + "outputs": [], + "source": [ + "pc = ParameterContainer()\n", + "rc = ReactionContainer()\n", + "pc.add([k4Vmax, k3r, k3f, k2f, j1pulse, PSOCE, J0_SOCE])\n", + "rc.add([r1, r2, r3, r4, rSOCE])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "15c35d39", + "metadata": {}, + "source": [ + "## Create/load in mesh\n", + "\n", + "In SMART we have different levels of meshes:\n", + "- Parent mesh: contains the entire geometry of the problem, including all surfaces and volumes\n", + "- Child meshes: submeshes (sections of the parent mesh) associated with individual compartments. Here, the child meshes are:\n", + " - Cyto: the portion of the outer cube outside of the inner cube, defined by `cell_markers = 1`\n", + " - ER: the inside portion of the inner cube, defined by `cell_markers = 2`\n", + " - PM: surface mesh where x=0, defined by `facet_markers = 10`\n", + " - ERm: surface mesh corresponding to all faces of the inner cube, defined by `facet_markers = 12`\n", + "\n", + "Here we create a UnitCube mesh as the Parent mesh, defined by\n", + "\n", + "$$\n", + "\\Omega = [0, 1] \\times [0, 1] \\times [0, 1] \\subset \\mathbb{R}^3\n", + "$$\n", + "\n", + "\n", + "The ER is a cube within the exterior cube, with dimensions 0.4 by 0.4 by 0.4 and a tunable gap between the ER and PM to test for different strengths of SOCE below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e77f97fe", + "metadata": {}, + "outputs": [], + "source": [ + "ER_PM_gap = 0 # must be greater than 0 and less than 0.6\n", + "ER_PM_gap = 0.05 * round(ER_PM_gap/0.05) # round to ensure flat edges in mesh\n", + "def cur_cube_condition(cell, xmin=ER_PM_gap, xmax=ER_PM_gap+0.4):\n", + " \"\"\"\n", + " Returns true when inside an inner cube region defined as:\n", + " xmin <= x <= xmax, ymin <= y <= ymax, zmin <= z <= zmax\n", + " \"\"\"\n", + " ymin = 0.3\n", + " ymax = 0.7\n", + " zmin = 0.3\n", + " zmax = 0.7\n", + " return (\n", + " (xmin - d.DOLFIN_EPS < cell.midpoint().x() < xmax + d.DOLFIN_EPS)\n", + " and (ymin - d.DOLFIN_EPS < cell.midpoint().y() < ymax + d.DOLFIN_EPS)\n", + " and (zmin - d.DOLFIN_EPS < cell.midpoint().z() < zmax + d.DOLFIN_EPS)\n", + " )\n", + "domain, facet_markers, cell_markers = mesh_tools.create_cubes(condition=cur_cube_condition, N=20)\n", + "visualization.plot_dolfin_mesh(domain, cell_markers, clip_plane=(1,\n", + " 1, 0), clip_origin=(0.5, 0.5, 0.5))\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "035113ed", + "metadata": {}, + "source": [ + "By default, `smart.mesh_tools.create_cubes` marks all faces of the outer cube as \"10\", our marker value associated with PM. Here, since we are only treating the x=0 face as PM, we alter the facet markers on all other faces, setting them equal to zero. They are then treated as no-flux boundaries not belonging to a designated surface compartment. The resultant mesh with the new facet and volume markers is displayed below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe56e162", + "metadata": {}, + "outputs": [], + "source": [ + "for face in d.faces(domain):\n", + " if face.midpoint().x() > d.DOLFIN_EPS and facet_markers[face] == 10:\n", + " facet_markers[face] = 0\n", + "img_mesh = mpimg.imread('example5-mesh.png')\n", + "plt.imshow(img_mesh)\n", + "plt.axis('off')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c17ebcbb", + "metadata": {}, + "source": [ + "We now save the mesh as an h5 file and then read it into SMART as a `ParentMesh` object. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5657aab1", + "metadata": {}, + "outputs": [], + "source": [ + "mesh_folder = pathlib.Path(\"mesh\")\n", + "mesh_folder.mkdir(exist_ok=True)\n", + "mesh_path = mesh_folder / \"DemoCuboidsMesh.h5\"\n", + "mesh_tools.write_mesh(\n", + " domain, facet_markers, cell_markers, filename=mesh_path\n", + ")\n", + "parent_mesh = mesh.ParentMesh(\n", + " mesh_filename=str(mesh_path),\n", + " mesh_filetype=\"hdf5\",\n", + " name=\"parent_mesh\",\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "0f1cf8a4", + "metadata": {}, + "source": [ + "## Model and solver initialization\n", + "\n", + "Now we are ready to set up the model. First we load the default configurations and set some configurations for the current solver." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c1a2a92", + "metadata": {}, + "outputs": [], + "source": [ + "conf = config.Config()\n", + "conf.solver.update(\n", + " {\n", + " \"final_t\": 1,\n", + " \"initial_dt\": 0.01,\n", + " \"time_precision\": 6,\n", + " }\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "1511acc7", + "metadata": {}, + "source": [ + "We create a model using the different containers and the parent mesh. For later reference, we save the model information as a pickle file. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3c3c30f", + "metadata": {}, + "outputs": [], + "source": [ + "model_cur = model.Model(pc, sc, cc, rc, conf, parent_mesh)\n", + "model_cur.to_pickle('model_cur.pkl')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "dfb70c74", + "metadata": {}, + "source": [ + "Note that we could later load the model information from the pickle file using the line:\n", + "```Python\n", + "model_cur = model.from_pickle(model_cur.pkl)\n", + "```\n", + "\n", + "Next we need to initialize the model and solver." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8976aa2", + "metadata": {}, + "outputs": [], + "source": [ + "model_cur.initialize()" + ] + }, + { + "cell_type": "markdown", + "id": "70f79003", + "metadata": {}, + "source": [ + "Now create mapping between PM and closest point on the ER surface for spatially dependent SOCE." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db8d9dca", + "metadata": {}, + "outputs": [], + "source": [ + "Pfun = PSOCE.dolfin_function # function associated with PSOCE values\n", + "V_PM = Pfun.function_space() # PM function space\n", + "# pull coordinate lists for PM and ERM\n", + "PMcoord = V_PM.tabulate_dof_coordinates()\n", + "ERMcoord = ERm.dolfin_mesh.coordinates()\n", + "# initialize vectors for distances and AER indices\n", + "AER_idx = np.zeros_like(Pfun.vector())\n", + "dist_vals = np.zeros_like(Pfun.vector())\n", + "# get parent mesh and vertex mapping to meshviews\n", + "mesh_ref = model_cur.parent_mesh.dolfin_mesh\n", + "ERM_map = ERm.dolfin_mesh.topology().mapping()[mesh_ref.id()].vertex_map()\n", + "ER_map = ER.dolfin_mesh.topology().mapping()[mesh_ref.id()].vertex_map()\n", + "# find closest ERM point for each PM coordinate\n", + "for i in range(len(PMcoord)):\n", + " curCoord = PMcoord[i]\n", + " dists = np.sqrt((curCoord[0]-ERMcoord[:,0])**2 + \n", + " (curCoord[1]-ERMcoord[:,1])**2 + \n", + " (curCoord[2]-ERMcoord[:,2])**2)\n", + " ERMidx = np.argmin(dists)\n", + " dist_vals[i] = dists[ERMidx]\n", + " global_idx = ERM_map[ERMidx]\n", + " ER_idx = np.nonzero(np.array(ER_map)==global_idx)[0]\n", + " if len(ER_idx) != 1:\n", + " raise ValueError(\"Node not found in ER mesh\")\n", + " AER_idx[i] = d.vertex_to_dof_map(AER.V)[ER_idx][0]\n", + "# now define PSOCE function values according to dist vals and AER values\n", + "d0 = 0.01\n", + "cref = 2.0\n", + "vals_new = (d0/(dist_vals+d0)) * (1/(1+(AER.u[\"u\"].vector()[AER_idx]/cref)**4))\n", + "PSOCE.dolfin_function.vector()[:] = vals_new\n", + "PSOCE.dolfin_function.vector().apply(\"insert\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d05a1a75", + "metadata": {}, + "source": [ + "## Solving system and storing data\n", + "\n", + "We create some XDMF files where we will store the output " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5df9509e", + "metadata": {}, + "outputs": [], + "source": [ + "# Write initial condition(s) to file\n", + "results = dict()\n", + "result_folder = pathlib.Path(f\"results_ER_PM_gap_{ER_PM_gap}\")\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", + "SOCEfile = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / f\"PSOCE.xdmf\"))\n", + "SOCEfile.parameters[\"flush_output\"] = True\n", + "SOCEfile.write(PSOCE.dolfin_function, model_cur.t)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "183376e6", + "metadata": {}, + "source": [ + "We now run the the solver and store the data at each time point to the initialized files. We also integrate A over the cytosolic volume at each time step to monitor the elevation in A over time and we display the concentration of A in the cytosol at t = 0.2 s." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31d42278", + "metadata": {}, + "outputs": [], + "source": [ + "# Set loglevel to warning in order not to pollute notebook output\n", + "logger.setLevel(logging.WARNING)\n", + "\n", + "avg_A = [A.initial_condition]\n", + "# define integration measure and total volume for computing average A at each time point\n", + "dx = d.Measure(\"dx\", domain=model_cur.cc['Cyto'].dolfin_mesh)\n", + "volume = d.assemble_mixed(1.0*dx)\n", + "# Solve\n", + "displayed = False\n", + "while model_cur.t < model_cur.final_t:\n", + " # Solve the system\n", + " model_cur.monolithic_solve()\n", + " # Update PSOCE\n", + " PSOCE.dolfin_function.vector()[:] = (d0/(dist_vals+d0)) * (1/(1+(AER.u[\"u\"].vector()[AER_idx]/cref)**4))\n", + " PSOCE.dolfin_function.vector().apply(\"insert\")\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", + " SOCEfile.write(PSOCE.dolfin_function, model_cur.t)\n", + " # compute average A concentration at each time step\n", + " int_val = d.assemble_mixed(model_cur.sc['A'].u['u']*dx)\n", + " avg_A.append(int_val / volume)\n", + " if model_cur.t >= 0.2 and not displayed:\n", + " visualization.plot(model_cur.sc['A'].u['u'],\n", + " clip_plane=(1, 1, 0), clip_origin=(0.5, 0.5, 0.5))\n", + " displayed = True" + ] + }, + { + "cell_type": "markdown", + "id": "621600d7", + "metadata": {}, + "source": [ + "Finally, we plot the average concentration of A over time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b8651e71", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(model_cur.tvec, avg_A)\n", + "plt.xlabel('Time (s)')\n", + "plt.ylabel('Cytosolic concentration of A (μM)')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d302bf3f", + "metadata": {}, + "source": [ + "We also compare the AUC for A with previous numerical simulations (regression test)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8b3c866", + "metadata": {}, + "outputs": [], + "source": [ + "tvec = np.zeros(len(model_cur.tvec))\n", + "for i in range(len(model_cur.tvec)):\n", + " tvec[i] = float(model_cur.tvec[i])\n", + "auc_cur = np.trapz(np.array(avg_A), tvec)\n", + "auc_compare = 4.646230684534995\n", + "percent_error = 100*np.abs(auc_cur - auc_compare)/auc_compare\n", + "assert percent_error < .01,\\\n", + " f\"Failed regression test: Example 5 results deviate {percent_error:.3f}% from the previous numerical solution\"" + ] + } + ], + "metadata": { + "jupytext": { + "cell_metadata_filter": "-all", + "main_language": "python", + "notebook_metadata_filter": "-all" + }, + "kernelspec": { + "display_name": "Python 3", + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/example5/model_cur.pkl b/examples/example5/model_cur.pkl new file mode 100644 index 0000000000000000000000000000000000000000..7830aa97204403756203f5688e83bd692ce48e07 GIT binary patch literal 3487 zcmcInO^h5z72cki{qg?U-mGGggCXFA$%+iNV~ZQvu(22K$m3k{@iFSBy< z{X$A+yq(3`^eX4m=H62fYa=R8kK+Ha zL-pvP`uL&x1SV&C>am`Byr({?o>ZL?2L1apXsd@Uhnf+b^gl8>?--r=KWQ?(doHhH zEv+`$r{;h0`X7FJx4GIhgZPs{6-$F~t0R(;WH#PMMBZ1+1NG>EH76~?QTcu+$O_p> zDgdB)6(hX5G_sLa$u=XDrmt{xT3Bqn{)#v8vU$nSi4n>?yNO+DZe}<6lF{2m7-HiY zjISBBNzH31M-Myxn%cU$`dz+DeMLQoy~)cpVA{}ZAsG1qV84cP-3uISFB(9%kU^vw zEZvLxTZT@K2!mI|prNVj?gavQ^|`qqOOmPt%x*|?fP=3v_pNS&rQr9$AJWYA=V`YEOvvz1Z*8M6C7`?h9W0?T4qS)`0MYNNp88TfF1xR2(V9IuQP6^i0yJ)(N zZqo^;`c@Uby8A0@m2ALpAl(oED>TcHm)e?u@uxb^8?`5D*iCL;zqWRXBgqBp#Khf^ zbRQfXj00|k780q|?R@K44nR1)an>I&T>bfhOl<~t_xF1=Qyg6jAeq3+- z#BTgdy~XybPp?K}r*Eguhr8Y#dGpQn)#aLa>_%!_a!#l(byNt8fkkCCenD$( zzs&B?>YL7X4}67s$7&^Me#B%0Pb*KoMT;S7Udmn&F<13&gHDdneu!(-l{?KqbAgOw zA7yM0NvS(;Rrb^`)o;2@_1oIgzq9)OzU!>)QDgI}0YT#~g4z1`JNG>9Rg#VsQaI?F zVK37t&LQi?iq)fL;hx^-j6mw*LJH|}PGsJOQI|B$)nQ*Rt0Q%!QBmXSmahZYBHaP8 zI_ha!7%=l)xtYpg<|2kCv#^XQuzmGqc{ts?5dCD?rt0QV` zoimbIQFH5DW(0F=?@oG%;m=Sr>!UcIt(Y~4)zR)eHr%ktZ#+~i4|aQYy=T{Y6@;4` z5$-{Rx?qJuQ3Ix|_T)X|1frJhVAYg0Lff*i12roU@1BAfK{W((Fta-0t>9PapRNK6XxEgn=D*HW*5O|_AX~L@-~h74mCE;8ufi^F((bP zw=Q#hbmh1wuJC~|H%WApqTtMD%R@hkrAvXQc>d{UWsZVa$7XRUA0b8 zM-J#TKe@%On_1}tPf{S%_jM@|(KI_KGGU35P^UFICc<4{%g17m7Hu~^(sfMH`0U++ Date: Mon, 23 Jun 2025 22:31:09 -0700 Subject: [PATCH 3/8] Remove model specs file for example 5 --- examples/example5/model_cur.pkl | Bin 3487 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 examples/example5/model_cur.pkl diff --git a/examples/example5/model_cur.pkl b/examples/example5/model_cur.pkl deleted file mode 100644 index 7830aa97204403756203f5688e83bd692ce48e07..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 3487 zcmcInO^h5z72cki{qg?U-mGGggCXFA$%+iNV~ZQvu(22K$m3k{@iFSBy< z{X$A+yq(3`^eX4m=H62fYa=R8kK+Ha zL-pvP`uL&x1SV&C>am`Byr({?o>ZL?2L1apXsd@Uhnf+b^gl8>?--r=KWQ?(doHhH zEv+`$r{;h0`X7FJx4GIhgZPs{6-$F~t0R(;WH#PMMBZ1+1NG>EH76~?QTcu+$O_p> zDgdB)6(hX5G_sLa$u=XDrmt{xT3Bqn{)#v8vU$nSi4n>?yNO+DZe}<6lF{2m7-HiY zjISBBNzH31M-Myxn%cU$`dz+DeMLQoy~)cpVA{}ZAsG1qV84cP-3uISFB(9%kU^vw zEZvLxTZT@K2!mI|prNVj?gavQ^|`qqOOmPt%x*|?fP=3v_pNS&rQr9$AJWYA=V`YEOvvz1Z*8M6C7`?h9W0?T4qS)`0MYNNp88TfF1xR2(V9IuQP6^i0yJ)(N zZqo^;`c@Uby8A0@m2ALpAl(oED>TcHm)e?u@uxb^8?`5D*iCL;zqWRXBgqBp#Khf^ zbRQfXj00|k780q|?R@K44nR1)an>I&T>bfhOl<~t_xF1=Qyg6jAeq3+- z#BTgdy~XybPp?K}r*Eguhr8Y#dGpQn)#aLa>_%!_a!#l(byNt8fkkCCenD$( zzs&B?>YL7X4}67s$7&^Me#B%0Pb*KoMT;S7Udmn&F<13&gHDdneu!(-l{?KqbAgOw zA7yM0NvS(;Rrb^`)o;2@_1oIgzq9)OzU!>)QDgI}0YT#~g4z1`JNG>9Rg#VsQaI?F zVK37t&LQi?iq)fL;hx^-j6mw*LJH|}PGsJOQI|B$)nQ*Rt0Q%!QBmXSmahZYBHaP8 zI_ha!7%=l)xtYpg<|2kCv#^XQuzmGqc{ts?5dCD?rt0QV` zoimbIQFH5DW(0F=?@oG%;m=Sr>!UcIt(Y~4)zR)eHr%ktZ#+~i4|aQYy=T{Y6@;4` z5$-{Rx?qJuQ3Ix|_T)X|1frJhVAYg0Lff*i12roU@1BAfK{W((Fta-0t>9PapRNK6XxEgn=D*HW*5O|_AX~L@-~h74mCE;8ufi^F((bP zw=Q#hbmh1wuJC~|H%WApqTtMD%R@hkrAvXQc>d{UWsZVa$7XRUA0b8 zM-J#TKe@%On_1}tPf{S%_jM@|(KI_KGGU35P^UFICc<4{%g17m7Hu~^(sfMH`0U++ Date: Mon, 23 Jun 2025 22:42:12 -0700 Subject: [PATCH 4/8] Change version number in pyproject --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 0bbb4e3..7a80519 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ requires = ["setuptools>=61.0.0", "wheel"] [project] name = "fenics-smart" -version = "2.2.3" +version = "2.3.0" description = "Spatial Modeling Algorithms for Reactions and Transport (SMART) is a high-performance finite-element-based simulation package for model specification and numerical simulation of spatially-varying reaction-transport processes in biological cells." authors = [{ name = "Justin Laughlin", email = "justinglaughlin@gmail.com" }] license = { file = "LICENSE" } From a25eabc6c17e1bdb5baf40bd802472e5d4f9ea70 Mon Sep 17 00:00:00 2001 From: emmetfrancis <99422170+emmetfrancis@users.noreply.github.com> Date: Tue, 24 Jun 2025 17:42:30 -0700 Subject: [PATCH 5/8] Update ER_PM_gap value and save Avals in example5 with SOCE --- examples/example5/example5_withSOCE.ipynb | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/example5/example5_withSOCE.ipynb b/examples/example5/example5_withSOCE.ipynb index d1405bd..c70075b 100644 --- a/examples/example5/example5_withSOCE.ipynb +++ b/examples/example5/example5_withSOCE.ipynb @@ -276,7 +276,7 @@ "A = Species(\"A\", 0.01, vol_unit, 1.0, D_unit, \"Cyto\")\n", "B = Species(\"B\", 0.0, vol_unit, 1.0, D_unit, \"Cyto\")\n", "AER = Species(\"AER\", 200.0, vol_unit, 5.0, D_unit, \"ER\")\n", - "# Uniform initial condition of R\n", + "# Uniform initial condition of R for simplicity\n", "R1 = Species(\"R1\", 1.0, surf_unit, 0.02, D_unit, \"ERm\")\n", "R1o = Species(\"R1o\", 0.0, surf_unit, 0.02, D_unit, \"ERm\")" ] @@ -490,7 +490,7 @@ "metadata": {}, "outputs": [], "source": [ - "ER_PM_gap = 0 # must be greater than 0 and less than 0.6\n", + "ER_PM_gap = 0.05 # must be greater than 0 and less than 0.6\n", "ER_PM_gap = 0.05 * round(ER_PM_gap/0.05) # round to ensure flat edges in mesh\n", "def cur_cube_condition(cell, xmin=ER_PM_gap, xmax=ER_PM_gap+0.4):\n", " \"\"\"\n", @@ -517,7 +517,7 @@ "id": "035113ed", "metadata": {}, "source": [ - "By default, `smart.mesh_tools.create_cubes` marks all faces of the outer cube as \"10\", our marker value associated with PM. Here, since we are only treating the x=0 face as PM, we alter the facet markers on all other faces, setting them equal to zero. They are then treated as no-flux boundaries not belonging to a designated surface compartment. The resultant mesh with the new facet and volume markers is displayed below." + "By default, `smart.mesh_tools.create_cubes` marks all faces of the outer cube as \"10\", our marker value associated with PM. Here, since we are only treating the x=0 face as PM, we alter the facet markers on all other faces, setting them equal to zero. They are then treated as no-flux boundaries not belonging to a designated surface compartment." ] }, { @@ -529,10 +529,7 @@ "source": [ "for face in d.faces(domain):\n", " if face.midpoint().x() > d.DOLFIN_EPS and facet_markers[face] == 10:\n", - " facet_markers[face] = 0\n", - "img_mesh = mpimg.imread('example5-mesh.png')\n", - "plt.imshow(img_mesh)\n", - "plt.axis('off')" + " facet_markers[face] = 0" ] }, { @@ -677,8 +674,8 @@ " raise ValueError(\"Node not found in ER mesh\")\n", " AER_idx[i] = d.vertex_to_dof_map(AER.V)[ER_idx][0]\n", "# now define PSOCE function values according to dist vals and AER values\n", - "d0 = 0.01\n", - "cref = 2.0\n", + "d0 = 0.05\n", + "cref = 100.0\n", "vals_new = (d0/(dist_vals+d0)) * (1/(1+(AER.u[\"u\"].vector()[AER_idx]/cref)**4))\n", "PSOCE.dolfin_function.vector()[:] = vals_new\n", "PSOCE.dolfin_function.vector().apply(\"insert\")" @@ -776,9 +773,13 @@ "metadata": {}, "outputs": [], "source": [ - "plt.plot(model_cur.tvec, avg_A)\n", + "plt.plot(model_cur.tvec, avg_A, label='PM-ER gap = 50 nm')\n", + "other_case = np.loadtxt('/root/shared/gitrepos/smart-dev/examples/example5/results_ER_PM_gap_0.2/Avals.txt')\n", + "plt.plot(other_case[0,:], other_case[1,:], label='PM-ER gap = 200 nm')\n", "plt.xlabel('Time (s)')\n", - "plt.ylabel('Cytosolic concentration of A (μM)')" + "plt.ylabel('Cytosolic concentration of A (μM)')\n", + "plt.legend()\n", + "np.savetxt(str(result_folder / f\"Avals.txt\"), [model_cur.tvec, avg_A])" ] }, { From 9ef3822fd50d7e4ea01cbfd36e16e3e58d942da7 Mon Sep 17 00:00:00 2001 From: emmetfrancis <99422170+emmetfrancis@users.noreply.github.com> Date: Sat, 12 Jul 2025 19:03:46 -0700 Subject: [PATCH 6/8] Add new meshing function for multicellular meshes --- smart/mesh_tools.py | 153 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 145 insertions(+), 8 deletions(-) diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py index fb04478..43456c1 100644 --- a/smart/mesh_tools.py +++ b/smart/mesh_tools.py @@ -719,6 +719,143 @@ def meshSizeCallback(dim, tag, x, y, z, lc): return (dmesh, mf2, mf3) +def create_multicell( + cubeSize: float = 100.0, + locVec: list = [[0, 0, 0]], + cellRad: float = 10.0, + hCube: float = 0, + hCell: float = 0, + interface_marker: int = 12, + outer_marker: int = 10, + extracell_tag: int = 2, + cell_vol_tag: int = 1, + comm: MPI.Comm = d.MPI.comm_world, + verbose: bool = False, +) -> Tuple[d.Mesh, d.MeshFunction, d.MeshFunction]: + """ + Creates a mesh with an outer cube containing embedded cells at specified locations. + Args: + cubeSize: Length of cube sides + locVec: vector of cell locations + hCube: maximum mesh size for cube + hCell: maximum mesh size for cell surfaces + interface_marker: The value to mark facets on the interface with + outer_marker: The value to mark facets on the outer ellipsoid with + inner_vol_tag: The value to mark the inner spherical volume with + outer_vol_tag: The value to mark the outer spherical volume with + comm: MPI communicator to create the mesh with + verbose: If true print gmsh output, else skip + Returns: + A triplet (mesh, facet_marker, cell_marker) + """ + import gmsh + + if np.isclose(cubeSize, 0): + ValueError("Outer cube size is equal to zero") + if np.isclose(hCube, 0): + hCube = 0.1 * max(cubeSize) + if np.isclose(hCell, 0): + hCell = 0.2 * cubeSize if np.isclose(cellRad, 0) else 0.2 * cellRad + # if innerRad > outerRad or innerLength >= outerLength: + # ValueError("Inner cylinder does not fit inside outer cylinder") + # Create the two cylinder mesh using gmsh + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", int(verbose)) + + gmsh.model.add("multicell") + # first add outer cube + cube = gmsh.model.occ.addBox( + -cubeSize / 2, -cubeSize / 2, -cubeSize / 2, cubeSize, cubeSize, cubeSize + ) + if np.isclose(cellRad, 0): + # Just a cube! + gmsh.model.occ.synchronize() + gmsh.model.add_physical_group(3, [cube], tag=extracell_tag) + facets = gmsh.model.getBoundary([(3, cube)]) + gmsh.model.add_physical_group(2, [facets[0][1]], tag=outer_marker) + else: + # Add cells + cell_list = [] + for i in range(len(locVec)): + cur_tag = gmsh.model.occ.addSphere(locVec[i][0], locVec[i][1], locVec[i][2], cellRad) + cell_list.append((3, cur_tag)) + # Create interface between cells and extracell + full_geo, maps = gmsh.model.occ.fragment([(3, cube)], cell_list) + cube_map = maps[0] + cell_maps = maps[1:] + gmsh.model.occ.synchronize() + + # Get the outer boundary + outer_shells = gmsh.model.getBoundary(full_geo, oriented=False) + # Get the inner boundary + inner_shells = [] + for i in range(len(cell_maps)): + inner_shells.append(gmsh.model.getBoundary(cell_maps[i], oriented=False)) + # Add physical markers for facets + gmsh.model.add_physical_group(2, [faces[1] for faces in outer_shells], tag=outer_marker) + gmsh.model.add_physical_group( + 2, [faces[0][1] for faces in inner_shells], tag=interface_marker + ) + + # Physical markers for + all_volumes = [tag[1] for tag in cube_map] + inner_volume = [tag[0][1] for tag in cell_maps] + outer_volume = [] + for vol in all_volumes: + if vol not in inner_volume: + outer_volume.append(vol) + gmsh.model.add_physical_group(3, outer_volume, tag=extracell_tag) + gmsh.model.add_physical_group(3, inner_volume, tag=cell_vol_tag) + + def meshSizeCallback(dim, tag, x, y, z, lc): + # mesh length is hEdge at the PM (defaults to 0.1*outerRad, + # or set when calling function) and hInnerEdge at the ERM + # (defaults to 0.2*innerRad, or set when calling function) + # between these, the value is interpolated based on r (polar coord), + # and inside the value is interpolated between hInnerEdge and 0.2*innerRad + # if innerRad=0, then the mesh length is interpolated between + # hEdge at the PM and 0.2*outerRad in the center + + if np.isclose(cellRad, 0): + return hCube + cell_locs = np.sqrt( + (x - np.array(locVec)[:, 0]) ** 2 + + (y - np.array(locVec)[:, 1]) ** 2 + + (z - np.array(locVec)[:, 2]) ** 2 + ) + closest_cell = min(cell_locs) + cellWeight = np.exp(-(closest_cell - cellRad) / (0.2 * cellRad)) + if closest_cell < cellRad: + return hCell + else: + return hCell * cellWeight + hCube * (1 - cellWeight) + + gmsh.model.mesh.setSizeCallback(meshSizeCallback) + # set off the other options for mesh size determination + gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + # this changes the algorithm from Frontal-Delaunay to Delaunay, + # which may provide better results when there are larger gradients in mesh size + gmsh.option.setNumber("Mesh.Algorithm", 5) + + gmsh.model.mesh.generate(3) + rank = MPI.COMM_WORLD.rank + tmp_folder = pathlib.Path(f"tmp_extracell_{cubeSize}_{cellRad}_{rank}") + tmp_folder.mkdir(exist_ok=True) + gmsh_file = tmp_folder / "extracell.msh" + gmsh.write(str(gmsh_file)) + gmsh.finalize() + + # return dolfin mesh of max dimension (parent mesh) and marker functions mf2 and mf3 + dmesh, mf2, mf3 = gmsh_to_dolfin(str(gmsh_file), tmp_folder, 3, comm) + # remove tmp mesh and tmp folder + gmsh_file.unlink(missing_ok=False) + tmp_folder.rmdir() + # return dolfin mesh, mf2 (2d tags) and mf3 (3d tags) + return (dmesh, mf2, mf3) + + def create_ellipses( xrad_outer: float = 3.0, yrad_outer: float = 1.0, @@ -1547,11 +1684,11 @@ def meshSizeCallback(dim, tag, x, y, z, lc): gmsh.finalize() # return dolfin mesh of max dimension (parent mesh) and marker functions mf2 and mf3 - dmesh, mf2, mf3 = gmsh_to_dolfin(str(gmsh_file), tmp_folder, 2, comm) + dmesh, mf1, mf2 = gmsh_to_dolfin(str(gmsh_file), tmp_folder, 2, comm) # remove tmp mesh and tmp folder gmsh_file.unlink(missing_ok=False) tmp_folder.rmdir() - # return dolfin mesh, mf2 (2d tags) and mf3 (3d tags) + # return dolfin mesh, mf1 (1d tags) and mf2 (2d tags) if return_curvature: if innerExpr == "": facet_list = [outer_marker] @@ -1560,7 +1697,7 @@ def meshSizeCallback(dim, tag, x, y, z, lc): facet_list = [outer_marker, interface_marker] cell_list = [outer_tag, inner_tag] if half_cell_with_curvature: # will likely not work in parallel... - dmesh_half, mf2_half, mf3_half = create_2Dcell( + dmesh_half, mf1_half, mf2_half = create_2Dcell( outerExpr, innerExpr, hEdge, @@ -1575,14 +1712,14 @@ def meshSizeCallback(dim, tag, x, y, z, lc): return_curvature=False, ) kappa_mf = compute_curvature( - dmesh, mf2, mf3, facet_list, cell_list, half_mesh_data=(dmesh_half, mf2_half) + dmesh, mf1, mf2, facet_list, cell_list, half_mesh_data=(dmesh_half, mf1_half) ) - (dmesh, mf2, mf3) = (dmesh_half, mf2_half, mf3_half) + (dmesh, mf1, mf2) = (dmesh_half, mf1_half, mf2_half) else: - kappa_mf = compute_curvature(dmesh, mf2, mf3, facet_list, cell_list) - return (dmesh, mf2, mf3, kappa_mf) + kappa_mf = compute_curvature(dmesh, mf1, mf2, facet_list, cell_list) + return (dmesh, mf1, mf2, kappa_mf) else: - return (dmesh, mf2, mf3) + return (dmesh, mf1, mf2) def gmsh_to_dolfin( From b1d24be32953cceb33ccb852aff7becdaad96bef Mon Sep 17 00:00:00 2001 From: emmetfrancis <99422170+emmetfrancis@users.noreply.github.com> Date: Mon, 14 Jul 2025 10:25:03 -0700 Subject: [PATCH 7/8] Update new example 7 for testing reaction-diffusion through multicellular network, ensure non-negativity of solution --- examples/example7/example7.ipynb | 304 +++++++++++++++++++++++++++++++ smart/mesh_tools.py | 101 +++++++--- 2 files changed, 379 insertions(+), 26 deletions(-) create mode 100644 examples/example7/example7.ipynb diff --git a/examples/example7/example7.ipynb b/examples/example7/example7.ipynb new file mode 100644 index 0000000..2c86caa --- /dev/null +++ b/examples/example7/example7.ipynb @@ -0,0 +1,304 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "f65f18d7", + "metadata": {}, + "source": [ + "# Example 7: Reaction-diffusion of molecule in a multicellular network\n", + "\n", + "Here, we consider the diffusion of a molecule from a central cell in a multicellular mesh to other cells in the environment." + ] + }, + { + "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, 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", + "\n", + "from matplotlib import pyplot as plt\n", + "import matplotlib.image as mpimg\n", + "from matplotlib import rcParams\n", + "\n", + "logger = logging.getLogger(\"smart\")\n", + "logger.setLevel(logging.INFO)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "95b9d865", + "metadata": {}, + "source": [ + "We define the relevant units here." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f4023cf", + "metadata": {}, + "outputs": [], + "source": [ + "# Aliases - base units\n", + "uM = unit.uM\n", + "um = unit.um\n", + "molecule = unit.molecule\n", + "sec = unit.sec\n", + "dimensionless = unit.dimensionless\n", + "# Aliases - units used in model\n", + "D_unit = um**2 / sec\n", + "flux_unit = uM * um / sec\n", + "vol_unit = uM\n", + "surf_unit = molecule / um**2" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "46582d26", + "metadata": {}, + "source": [ + "## Model generation\n", + "\n", + "We define the compartments and species first, with their respective containers." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02a000f2", + "metadata": {}, + "outputs": [], + "source": [ + "EC = Compartment(\"EC\", 3, um, 1)\n", + "PM1 = Compartment(\"PM1\", 2, um, 11) # source cell\n", + "PM2 = Compartment(\"PM2\", 2, um, 12) # other cells\n", + "PM1.specify_nonadjacency(['PM2'])\n", + "PM2.specify_nonadjacency(['PM1'])\n", + "\n", + "cc = CompartmentContainer()\n", + "cc.add([EC, PM1, PM2])\n", + "\n", + "G = Species(\"G\", 1.0, vol_unit, 1000.0, D_unit, \"EC\")\n", + "Rbound = Species(\"Rbound\", 0.0, surf_unit, 0.0, D_unit, \"PM2\")\n", + "sc = SpeciesContainer()\n", + "sc.add([G, Rbound])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "3c56e840", + "metadata": {}, + "source": [ + "Define parameters and reactions, then place in respective containers. Here, there are 3 reactions:\n", + "* r1: release of G from PM1\n", + "* r2: binding of G to PM2\n", + "* r3: degradation of G" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e1f6882", + "metadata": {}, + "outputs": [], + "source": [ + "# G release from cell 1\n", + "j1pulse = Parameter(\"j1pulse\", 1000.0, flux_unit)\n", + "r1 = Reaction(\n", + " \"r1\",\n", + " [],\n", + " [\"G\"],\n", + " param_map={\"J\": \"j1pulse\"},\n", + " eqn_f_str=\"J\",\n", + " explicit_restriction_to_domain=\"PM1\",\n", + ")\n", + "# G binding to PM2\n", + "Rtot = Parameter(\"Rtot\", 100.0, surf_unit)\n", + "kbind = Parameter(\"kbind\", 1.0, 1/(uM*sec))\n", + "kunbind = Parameter(\"kunbind\", 0.01, 1/sec)\n", + "r2 = Reaction(\"r2\", [\"G\"], [\"Rbound\"],\n", + " param_map={\"on\":\"kbind\",\"off\":\"kunbind\",\"Rtot\":\"Rtot\"},\n", + " eqn_f_str=\"on*G*(Rtot-Rbound) - off*Rbound\",\n", + " explicit_restriction_to_domain=\"PM2\")\n", + "\n", + "# G degradation\n", + "kdeg = Parameter(\"kdeg\", 0.01, 1/sec)\n", + "r3 = Reaction(\"r3\", [\"G\"], [], param_map={\"k\":\"kdeg\"},\n", + " eqn_f_str=\"k*G\")\n", + "\n", + "pc = ParameterContainer()\n", + "pc.add([j1pulse,Rtot,kbind,kunbind,kdeg])\n", + "rc = ReactionContainer()\n", + "rc.add([r1,r2,r3])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "15c35d39", + "metadata": {}, + "source": [ + "## Create and load in mesh\n", + "\n", + "Here, we consider cells embedded in a cube mesh. The source cell is located at (0,0,0) and 8 other cells are spread equidistant through the mesh." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe56e162", + "metadata": {}, + "outputs": [], + "source": [ + "domain, facet_markers, cell_markers = mesh_tools.create_multicell(cubeSize=100.0, locVec1=[[0,0,0]],\n", + " locVec2=[[-30,-30,-30], [-30,30,-30], [-30,-30,30], [-30,30,30], \n", + " [30,-30,-30], [30,30,-30], [30,-30,30], [30,30,30]],\n", + " cellRad1 = 10.0, cellRad2 = 10.0, hCube = 5.0, hCell = 2.0)\n", + "# Write mesh and meshfunctions to file\n", + "mesh_folder = pathlib.Path(\"mesh\")\n", + "mesh_folder.mkdir(exist_ok=True)\n", + "mesh_path = mesh_folder / \"multicell_mesh.h5\"\n", + "mesh_tools.write_mesh(\n", + " domain, facet_markers, cell_markers, filename=mesh_path\n", + ")\n", + "parent_mesh = mesh.ParentMesh(\n", + " mesh_filename=str(mesh_path),\n", + " mesh_filetype=\"hdf5\",\n", + " name=\"parent_mesh\",\n", + ")\n", + "visualization.plot_dolfin_mesh(domain, cell_markers, facet_markers)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "0943588e", + "metadata": {}, + "source": [ + "Initialize model and solver." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac88bdec", + "metadata": {}, + "outputs": [], + "source": [ + "config_cur = config.Config()\n", + "config_cur.flags.update({\"allow_unused_components\": True})\n", + "model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh)\n", + "config_cur.solver.update(\n", + " {\n", + " \"final_t\": 10.0,\n", + " \"initial_dt\": 0.001,\n", + " \"time_precision\": 8,\n", + " \"reset_timestep_for_negative_solution\": True,\n", + " }\n", + ")\n", + "model_cur.initialize()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "5d5aacbd", + "metadata": {}, + "source": [ + "Initialize XDMF files for saving results, save model information to .pkl file, then solve the system until `model_cur.t > model_cur.final_t`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b54d28ca", + "metadata": {}, + "outputs": [], + "source": [ + "# Write initial condition(s) to file\n", + "results = dict()\n", + "result_folder = pathlib.Path(f\"results\")\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", + "model_cur.to_pickle(\"model_cur.pkl\")\n", + "\n", + "# Set loglevel to warning in order not to pollute notebook output\n", + "logger.setLevel(logging.WARNING)\n", + "# Solve\n", + "displayed = False\n", + "while True:\n", + " # Solve the system\n", + " model_cur.monolithic_solve()\n", + " model_cur.adjust_dt()\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", + " print(f\"Done with t={model_cur.t}\")\n", + " # End if we've passed the final time\n", + " if model_cur.t >= model_cur.final_t:\n", + " break" + ] + } + ], + "metadata": { + "jupytext": { + "cell_metadata_filter": "-all", + "main_language": "python", + "notebook_metadata_filter": "-all" + }, + "kernelspec": { + "display_name": "Python 3", + "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 +} diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py index 43456c1..4ffabef 100644 --- a/smart/mesh_tools.py +++ b/smart/mesh_tools.py @@ -721,14 +721,18 @@ def meshSizeCallback(dim, tag, x, y, z, lc): def create_multicell( cubeSize: float = 100.0, - locVec: list = [[0, 0, 0]], - cellRad: float = 10.0, + locVec1: list = [[0, 0, 0]], + locVec2: list = [], + cellRad1: float = 10.0, + cellRad2: float = 10.0, hCube: float = 0, hCell: float = 0, - interface_marker: int = 12, + interface_marker1: int = 11, + interface_marker2: int = 12, outer_marker: int = 10, - extracell_tag: int = 2, - cell_vol_tag: int = 1, + extracell_tag: int = 1, + cell_vol_tag1: int = 2, + cell_vol_tag2: int = 3, comm: MPI.Comm = d.MPI.comm_world, verbose: bool = False, ) -> Tuple[d.Mesh, d.MeshFunction, d.MeshFunction]: @@ -755,7 +759,7 @@ def create_multicell( if np.isclose(hCube, 0): hCube = 0.1 * max(cubeSize) if np.isclose(hCell, 0): - hCell = 0.2 * cubeSize if np.isclose(cellRad, 0) else 0.2 * cellRad + hCell = 0.2 * cubeSize if np.isclose(cellRad1, 0) else 0.2 * cellRad1 # if innerRad > outerRad or innerLength >= outerLength: # ValueError("Inner cylinder does not fit inside outer cylinder") # Create the two cylinder mesh using gmsh @@ -767,7 +771,9 @@ def create_multicell( cube = gmsh.model.occ.addBox( -cubeSize / 2, -cubeSize / 2, -cubeSize / 2, cubeSize, cubeSize, cubeSize ) - if np.isclose(cellRad, 0): + if (np.isclose(cellRad1, 0) or len(locVec1) == 0) and ( + np.isclose(cellRad2, 0) or len(locVec2) == 0 + ): # Just a cube! gmsh.model.occ.synchronize() gmsh.model.add_physical_group(3, [cube], tag=extracell_tag) @@ -776,8 +782,17 @@ def create_multicell( else: # Add cells cell_list = [] - for i in range(len(locVec)): - cur_tag = gmsh.model.occ.addSphere(locVec[i][0], locVec[i][1], locVec[i][2], cellRad) + # first add source cell(s) + for i in range(len(locVec1)): + cur_tag = gmsh.model.occ.addSphere( + locVec1[i][0], locVec1[i][1], locVec1[i][2], cellRad1 + ) + cell_list.append((3, cur_tag)) + # now add additional cells + for i in range(len(locVec2)): + cur_tag = gmsh.model.occ.addSphere( + locVec2[i][0], locVec2[i][1], locVec2[i][2], cellRad2 + ) cell_list.append((3, cur_tag)) # Create interface between cells and extracell full_geo, maps = gmsh.model.occ.fragment([(3, cube)], cell_list) @@ -788,24 +803,32 @@ def create_multicell( # Get the outer boundary outer_shells = gmsh.model.getBoundary(full_geo, oriented=False) # Get the inner boundary - inner_shells = [] - for i in range(len(cell_maps)): - inner_shells.append(gmsh.model.getBoundary(cell_maps[i], oriented=False)) + inner_shells1 = [] + inner_shells2 = [] + for i in range(0, len(locVec1)): + inner_shells1.append(gmsh.model.getBoundary(cell_maps[i], oriented=False)) + for i in range(len(locVec1), len(cell_maps)): + inner_shells2.append(gmsh.model.getBoundary(cell_maps[i], oriented=False)) # Add physical markers for facets gmsh.model.add_physical_group(2, [faces[1] for faces in outer_shells], tag=outer_marker) gmsh.model.add_physical_group( - 2, [faces[0][1] for faces in inner_shells], tag=interface_marker + 2, [faces[0][1] for faces in inner_shells1], tag=interface_marker1 + ) + gmsh.model.add_physical_group( + 2, [faces[0][1] for faces in inner_shells2], tag=interface_marker2 ) # Physical markers for all_volumes = [tag[1] for tag in cube_map] - inner_volume = [tag[0][1] for tag in cell_maps] + inner_volumes1 = [tag[0][1] for tag in cell_maps[0 : len(locVec1)]] + inner_volumes2 = [tag[0][1] for tag in cell_maps[len(locVec1) :]] outer_volume = [] for vol in all_volumes: - if vol not in inner_volume: + if (vol not in inner_volumes1) and (vol not in inner_volumes2): outer_volume.append(vol) gmsh.model.add_physical_group(3, outer_volume, tag=extracell_tag) - gmsh.model.add_physical_group(3, inner_volume, tag=cell_vol_tag) + gmsh.model.add_physical_group(3, inner_volumes1, tag=cell_vol_tag1) + gmsh.model.add_physical_group(3, inner_volumes2, tag=cell_vol_tag2) def meshSizeCallback(dim, tag, x, y, z, lc): # mesh length is hEdge at the PM (defaults to 0.1*outerRad, @@ -816,18 +839,44 @@ def meshSizeCallback(dim, tag, x, y, z, lc): # if innerRad=0, then the mesh length is interpolated between # hEdge at the PM and 0.2*outerRad in the center - if np.isclose(cellRad, 0): + if (np.isclose(cellRad1, 0) or len(locVec1) == 0) and ( + np.isclose(cellRad2, 0) or len(locVec2) == 0 + ): return hCube - cell_locs = np.sqrt( - (x - np.array(locVec)[:, 0]) ** 2 - + (y - np.array(locVec)[:, 1]) ** 2 - + (z - np.array(locVec)[:, 2]) ** 2 - ) - closest_cell = min(cell_locs) - cellWeight = np.exp(-(closest_cell - cellRad) / (0.2 * cellRad)) - if closest_cell < cellRad: + elif np.isclose(cellRad1, 0) or len(locVec1) == 0: + cell_locs1 = [np.inf] + cell_locs2 = np.sqrt( + (x - np.array(locVec2)[:, 0]) ** 2 + + (y - np.array(locVec2)[:, 1]) ** 2 + + (z - np.array(locVec2)[:, 2]) ** 2 + ) + elif np.isclose(cellRad2, 0) or len(locVec2) == 0: + cell_locs1 = np.sqrt( + (x - np.array(locVec1)[:, 0]) ** 2 + + (y - np.array(locVec1)[:, 1]) ** 2 + + (z - np.array(locVec1)[:, 2]) ** 2 + ) + cell_locs2 = [np.inf] + else: + cell_locs1 = np.sqrt( + (x - np.array(locVec1)[:, 0]) ** 2 + + (y - np.array(locVec1)[:, 1]) ** 2 + + (z - np.array(locVec1)[:, 2]) ** 2 + ) + cell_locs2 = np.sqrt( + (x - np.array(locVec2)[:, 0]) ** 2 + + (y - np.array(locVec2)[:, 1]) ** 2 + + (z - np.array(locVec2)[:, 2]) ** 2 + ) + closest_cell1 = min(cell_locs1) + closest_cell2 = min(cell_locs2) + if (closest_cell1 < cellRad1) or (closest_cell2 < cellRad2): return hCell else: + if closest_cell1 < closest_cell2: + cellWeight = np.exp(-(closest_cell1 - cellRad1) / (0.2 * cellRad1)) + else: + cellWeight = np.exp(-(closest_cell2 - cellRad2) / (0.2 * cellRad2)) return hCell * cellWeight + hCube * (1 - cellWeight) gmsh.model.mesh.setSizeCallback(meshSizeCallback) @@ -841,7 +890,7 @@ def meshSizeCallback(dim, tag, x, y, z, lc): gmsh.model.mesh.generate(3) rank = MPI.COMM_WORLD.rank - tmp_folder = pathlib.Path(f"tmp_extracell_{cubeSize}_{cellRad}_{rank}") + tmp_folder = pathlib.Path(f"tmp_extracell_{cubeSize}_{cellRad1}_{cellRad2}_{rank}") tmp_folder.mkdir(exist_ok=True) gmsh_file = tmp_folder / "extracell.msh" gmsh.write(str(gmsh_file)) From 25e70770bc1f1b98fa846ec0a2f6acd397ad37bc Mon Sep 17 00:00:00 2001 From: emmetfrancis <99422170+emmetfrancis@users.noreply.github.com> Date: Mon, 14 Jul 2025 10:32:20 -0700 Subject: [PATCH 8/8] bump version number in pyproject --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 7a80519..e444256 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ requires = ["setuptools>=61.0.0", "wheel"] [project] name = "fenics-smart" -version = "2.3.0" +version = "2.3.0.beta.1" description = "Spatial Modeling Algorithms for Reactions and Transport (SMART) is a high-performance finite-element-based simulation package for model specification and numerical simulation of spatially-varying reaction-transport processes in biological cells." authors = [{ name = "Justin Laughlin", email = "justinglaughlin@gmail.com" }] license = { file = "LICENSE" }