Example 2 - 3D: Simple 3D cell signaling model#
We model a reaction between the cell interior and cell membrane within a dendritic spine:
Cyto - 3D spine volume
PM - 2D cell boundary
Model from Rangamani et al, 2013, Cell. A cytosolic species, “A”, reacts with a species on the PM, “B”, to form a new species on the PM, “X”.
\[\begin{split}
\frac{\partial{C_A}}{\partial{t}} = D_A \nabla ^2 C_A \quad \text{in} \; \Omega_{Cyto}\\
\text{B.C. on A:} \quad D_A (\textbf{n} \cdot \nabla C_A) = -k_{on} C_A N_X + k_{off} N_B \quad \text{on} \; \Gamma_{PM}
\end{split}\]
Similarly, the PDEs for X and B are given by: $\( \frac{\partial{N_X}}{\partial{t}} = D_X \nabla ^2 N_X - k_{on} C_A N_X + k_{off} N_B \quad \text{on} \; \Gamma_{PM}\\ \frac{\partial{N_B}}{\partial{t}} = D_B \nabla ^2 N_B + k_{on} C_A N_X - k_{off} N_B \quad \text{on} \; \Gamma_{PM} \)$
import matplotlib.image as mpimg
from matplotlib import pyplot as plt
img_A = mpimg.imread('axb-diagram.png')
plt.imshow(img_A)
plt.axis('off')
(-0.5, 5830.5, 3140.5, -0.5)
Imports and logger initialization:
import logging
import pathlib
import dolfin as d
import numpy as np
from smart import config, mesh, mesh_tools, model, visualization
from smart.model_assembly import (Compartment, CompartmentContainer, Parameter,
ParameterContainer, Reaction,
ReactionContainer, Species, SpeciesContainer)
from smart.units import unit
logger = logging.getLogger("smart")
logger.setLevel(logging.INFO)
/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4
warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
First, we define the various units for use in the model.
um = unit.um
molecule = unit.molecule
sec = unit.sec
dimensionless = unit.dimensionless
D_unit = um**2 / sec
vol_unit = unit.uM
flux_unit = molecule / (um**2 * sec)
surf_unit = molecule / um**2
Next we generate the model by assembling the compartment, species, parameter, and reaction containers (see Example 1 for or API documentation for more details).
# =============================================================================================
# Compartments
# =============================================================================================
# name, topological dimensionality, length scale units, marker value
Cyto = Compartment("Cyto", 3, um, 1)
PM = Compartment("PM", 2, um, 10)
cc = CompartmentContainer()
cc.add([Cyto, PM])
# =============================================================================================
# Species
# =============================================================================================
# name, initial concentration, concentration units, diffusion, diffusion units, compartment
A = Species("A", 1.0, vol_unit, 1.0, D_unit, "Cyto")
X = Species("X", 1000, surf_unit, 0.1, D_unit, "PM")
B = Species("B", 0.0, surf_unit, 0.01, D_unit, "PM")
sc = SpeciesContainer()
sc.add([A, X, B])
# =============================================================================================
# Parameters and Reactions
# =============================================================================================
# Reaction of A and X to make B (Cyto-PM reaction)
kon = Parameter("kon", 1.0, 1/(vol_unit*sec))
koff = Parameter("koff", 0.1, 1/sec)
r1 = Reaction("r1", ["A", "X"], ["B"],
param_map={"on": "kon", "off": "koff"},
species_map={"A": "A", "X": "X", "B": "B"})
pc = ParameterContainer()
pc.add([kon, koff])
rc = ReactionContainer()
rc.add([r1])
Now we load in the dendritic spine mesh and generate the marker functions mf3 and mf2.
# Load mesh
spine_mesh = d.Mesh('spine_mesh.xml')
mf3 = d.MeshFunction("size_t", spine_mesh, 3, 1)
mf2 = d.MeshFunction("size_t", spine_mesh, 2, spine_mesh.domains())
visualization.plot_dolfin_mesh(spine_mesh, mf3, clip_origin=(0.143, 0.107, -0.065))
Write mesh and meshfunctions to file, then create mesh.ParentMesh object.
mesh_folder = pathlib.Path("spine_mesh")
mesh_folder.mkdir(exist_ok=True)
mesh_file = mesh_folder / "spine_mesh.h5"
mesh_tools.write_mesh(spine_mesh, mf2, mf3, mesh_file)
parent_mesh = mesh.ParentMesh(
mesh_filename=str(mesh_file),
mesh_filetype="hdf5",
name="parent_mesh",
)
2024-05-07 12:02:08,100 smart.mesh - INFO - HDF5 mesh, "parent_mesh", successfully loaded from file: spine_mesh/spine_mesh.h5! (mesh.py:245)
2024-05-07 12:02:08,101 smart.mesh - INFO - 0 subdomains successfully loaded from file: spine_mesh/spine_mesh.h5! (mesh.py:258)
Initialize model and solvers.
configCur = config.Config()
configCur.solver.update(
{
"final_t": 5.0,
"initial_dt": 0.05,
"time_precision": 6,
"use_snes": True,
}
)
modelCur = model.Model(pc, sc, cc, rc, configCur, parent_mesh)
modelCur.initialize()
2024-05-07 12:02:19,226 smart.model - WARNING - Warning! Initial L2-norm of compartment PM is 111.62097871138167 (possibly too large). (model.py:1179)
2024-05-07 12:02:19,618 smart.solvers - INFO - Jpetsc_nest assembled, size = (24802, 24802) (solvers.py:199)
2024-05-07 12:02:19,619 smart.solvers - INFO - Initializing block residual vector (solvers.py:207)
2024-05-07 12:02:19,908 smart.model_assembly - INFO -
╒══════╤═════════════════╤═════════════════════╤════════════╤═════════╤═════════╕
│ │ quantity │ is_time_dependent │ sym_expr │ notes │ group │
╞══════╪═════════════════╪═════════════════════╪════════════╪═════════╪═════════╡
│ kon │ 1.00×10⁰ 1/s/µM │ False │ │ │ │
├──────┼─────────────────┼─────────────────────┼────────────┼─────────┼─────────┤
│ koff │ 1.00×10⁻¹ 1/s │ False │ │ │ │
╘══════╧═════════════════╧═════════════════════╧════════════╧═════════╧═════════╛
(model_assembly.py:337)
2024-05-07 12:02:19,913 smart.model_assembly - INFO -
╒════╤════════════════════╤═════════════════╤═════════════════════╤═══════════════════════╕
│ │ compartment_name │ Diffusion │ initial_condition │ concentration_units │
╞════╪════════════════════╪═════════════════╪═════════════════════╪═══════════════════════╡
│ A │ Cyto │ 1.00×10⁰ µm²/s │ 1 │ 1.00×10⁰ µM │
├────┼────────────────────┼─────────────────┼─────────────────────┼───────────────────────┤
│ X │ PM │ 1.00×10⁻¹ µm²/s │ 1000 │ 1.00×10⁰ molecule/µm² │
├────┼────────────────────┼─────────────────┼─────────────────────┼───────────────────────┤
│ B │ PM │ 1.00×10⁻² µm²/s │ 0 │ 1.00×10⁰ molecule/µm² │
╘════╧════════════════════╧═════════════════╧═════════════════════╧═══════════════════════╛
(model_assembly.py:337)
2024-05-07 12:02:20,667 smart.model_assembly - INFO -
╒══════╤══════════════════╤═══════════════╤════════════════╤════════════╤═════════════╤═══════════════╤═══════════════╕
│ │ dimensionality │ num_species │ num_vertices │ num_dofs │ num_cells │ cell_marker │ nvolume │
╞══════╪══════════════════╪═══════════════╪════════════════╪════════════╪═════════════╪═══════════════╪═══════════════╡
│ Cyto │ 3 │ 1 │ 11562 │ 11562 │ 48096 │ 1 │ 6.71×10⁻¹ µm³ │
├──────┼──────────────────┼───────────────┼────────────────┼────────────┼─────────────┼───────────────┼───────────────┤
│ PM │ 2 │ 2 │ 6620 │ 13240 │ 13107 │ 10 │ 6.24×10⁰ µm² │
╘══════╧══════════════════╧═══════════════╧════════════════╧════════════╧═════════════╧═══════════════╧═══════════════╛
(model_assembly.py:337)
2024-05-07 12:02:20,672 smart.model_assembly - INFO -
╒════╤════════════╤═══════╤═════════════╤═════════════╕
│ │ lhs │ rhs │ eqn_f_str │ eqn_r_str │
╞════╪════════════╪═══════╪═════════════╪═════════════╡
│ r1 │ ['A', 'X'] │ ['B'] │ kon*A*X │ koff*B │
╘════╧════════════╧═══════╧═════════════╧═════════════╛
(model_assembly.py:337)
╒════╤═════════════╤══════╤══════════════════╤═════════════╤══════════════╤════════════════╕
│ │ name │ id │ dimensionality │ num_cells │ num_facets │ num_vertices │
╞════╪═════════════╪══════╪══════════════════╪═════════════╪══════════════╪════════════════╡
│ 0 │ parent_mesh │ 10 │ 3 │ 48096 │ 103186 │ 11562 │
├────┼─────────────┼──────┼──────────────────┼─────────────┼──────────────┼────────────────┤
│ 1 │ Cyto │ 24 │ 3 │ 48096 │ 103186 │ 11562 │
├────┼─────────────┼──────┼──────────────────┼─────────────┼──────────────┼────────────────┤
│ 2 │ PM │ 31 │ 2 │ 13107 │ 19727 │ 6620 │
╘════╧═════════════╧══════╧══════════════════╧═════════════╧══════════════╧════════════════╛
Save model information to .pkl file and write initial conditions to file.
modelCur.to_pickle('modelCur.pkl')
results = dict()
result_folder = pathlib.Path("resultsSpine")
result_folder.mkdir(exist_ok=True)
for species_name, species in modelCur.sc.items:
results[species_name] = d.XDMFFile(
modelCur.mpi_comm_world, str(result_folder / f"{species_name}.xdmf")
)
results[species_name].parameters["flush_output"] = True
results[species_name].write(modelCur.sc[species_name].u["u"], modelCur.t)
Solve the system until modelCur.t > modelCur.final_t. We display the surface distribution of species B at \(t\) = 1.0 s for comparison with Figure 1 of the JOSS paper.
tvec = [0]
avg_B = [B.initial_condition]
# Set loglevel to warning in order not to pollute notebook output
logger.setLevel(logging.WARNING)
displayed = False
while True:
# Solve the system
modelCur.monolithic_solve()
# Save results for post processing
for species_name, species in modelCur.sc.items:
results[species_name].write(modelCur.sc[species_name].u["u"], modelCur.t)
dx = d.Measure("dx", domain=modelCur.cc['PM'].dolfin_mesh)
int_val = d.assemble_mixed(modelCur.sc['B'].u['u']*dx)
volume = d.assemble_mixed(1.0*dx)
avg_B.append(int_val / volume)
tvec.append(modelCur.t)
if modelCur.t >= 1.0 and not displayed:
visualization.plot(modelCur.sc["B"].u["u"], clip_logic=False, clim=(50,75))
displayed = True
# End if we've passed the final time
if modelCur.t >= modelCur.final_t:
break
Now we plot the average concentration of B in the dendritic spine over time, which should match the plot shown in Fig 1 of JOSS.
avg_B = np.array(avg_B)
plt.plot(tvec, 1e12*avg_B*1e12/(6.02e23), label='SMART simulation')
plt.xlabel('Time (s)')
plt.ylabel('B $\\times$ $\mathrm{10^{12}~(mol/m^2)}$')
Text(0, 0.5, 'B $\\times$ $\\mathrm{10^{12}~(mol/m^2)}$')
Regression test against previous runs.
percent_error = 100*np.abs(avg_B[-1]-64.35444884201829)
assert percent_error < .01,\
f"Failed regression test: Example 2-3d results deviate {percent_error:.3f}% from the previous numerical solution"