# Spatial modeling algorithms for reactions and transport (SMART) in
biological cells

Emmet A. Francis  Department of Mechanical and Aerospace Engineering,
University of California San Diego, La Jolla, CA, USA  these authors
contributed equally to this work  Justin G. Laughlin  Department of Mechanical
and Aerospace Engineering, University of California San Diego, La Jolla, CA,
USA  Computational Engineering Division, Lawrence Livermore National
Laboratory, Livermore, CA, USA  these authors contributed equally to this work
Jørgen S. Dokken  Department of Numerical Analysis and Scientific Computing,
Simula Research Laboratory, Oslo, Norway  Henrik N. T. Finsberg  Department of
Computational Physiology, Simula Research Laboratory, Oslo, Norway
Christopher T. Lee  Department of Mechanical and Aerospace Engineering,
University of California San Diego, La Jolla, CA, USA  Marie E. Rognes
Department of Numerical Analysis and Scientific Computing, Simula Research
Laboratory, Oslo, Norway<EMAIL_ADDRESS><EMAIL_ADDRESS>Padmini
Rangamani  Department of Mechanical and Aerospace Engineering, University of
California San Diego, La Jolla, CA, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Biological cells rely on precise spatiotemporal coordination of biochemical
reactions to control their many functions. Such cell signaling networks have
been a common focus for mathematical models, but they remain challenging to
simulate, particularly in realistic cell geometries. Herein, we present our
software, Spatial Modeling Algorithms for Reaction and Transport (SMART), a
package that takes in high-level user specifications about cell signaling
networks and molecular transport, and then assembles and solves the associated
mathematical and computational systems. SMART uses state-of-the-art finite
element analysis, via the FEniCS Project software, to efficiently and
accurately resolve cell signaling events over discretized cellular and
subcellular geometries. We demonstrate its application to several different
biological systems, including YAP/TAZ mechanotransduction, calcium signaling
in neurons and cardiomyocytes, and ATP generation in mitochondria. Throughout,
we utilize experimentally-derived realistic cellular geometries represented by
well-conditioned tetrahedral meshes. These scenarios demonstrate the
applicability, flexibility, accuracy and efficiency of SMART across a range of
temporal and spatial scales.

## Introduction

In the past decade, computational modeling has become an integral part of the
discovery toolkit in biology along with advances in experimental technologies
[1, 2, 3]. One of the fundamental tenets of biology is that structure and
function are very closely related [4, 5]. In single cell biology, this is
reflected by the spatial compartmentalization of cellular signaling in
different subcellular locations and organelles. Advances in microscopy in
recent decades have provided data to support this notion from two approaches:
electron microscopy for a detailed characterization of subcellular structure
[6, 7, 8, 9] and super-resolution microscopy for the spatiotemporal
localization of biochemical species that are important for cellular functions
including signaling [10, 11]. Detailed computational models using realistic
cellular geometries and reaction-diffusion mechanisms can help us identify the
possible biophysical mechanisms underlying such spatial compartmentalization
[12, 13, 14, 15]. However, representing these details including subcellular
geometries such as organelles and the relevant reaction-transport formulations
using the appropriate computational description remains an open challenge.

Historically, many mathematical models of cell signaling have neglected
spatial effects, treating the cell as a well-mixed volume. In certain cases,
this approximation can be justified, but given the slow diffusion of certain
signaling molecules, the crowded intracellular environment, and the complexity
of cellular geometries, such approximations can diminish the predictive
capability of the models. Furthermore, due to the variety of membrane-bound
organelles present in cells, reaction networks are coupled across subvolumes
and involve a complicated set of bulk-surface reactions at organelle membranes
(Figure 1A). As a result, moving from well-mixed models to multicompartment
spatiotemporal models of cell signaling presents many technical barriers.
Mathematically, this involves transitioning from systems of ordinary
differential equations (ODEs) to systems of mixed-dimensional partial
differential equations (PDEs). Such PDE systems are notoriously difficult to
solve numerically due to nonlinearities, stiffness, and instabilities [16].
Furthermore, solving these equations within realistic cellular geometries
requires robust discretization of complex geometries [8, 9] (Figure 1B-C) and
is computationally expensive due to the high dimensionality of the systems.

Recent efforts have been made to accurately represent cellular and subcellular
geometries with well-conditioned triangular and tetrahedral meshes [17, 18,
19, 20]. In particular, GAMer2 (Geometry-preserving Adaptive MeshER version 2)
allows users to convert microscopy images of cells into suitable meshes for
finite element simulations [19, 20]. These meshes can be annotated to mark the
location of subcellular structures such as the nucleus, endoplasmic reticulum
(ER), or mitochondria, along with their respective membrane boundaries [21,
22]. Meshes from GAMer2 can be readily loaded into the open-source finite
element software package, FEniCS [23], and used to simulate spatial models of
cell signaling in realistic geometries. However, translating systems of PDEs
defining cell signaling systems into variational forms defined over
subsurfaces and subvolumes of the geometry, and solving these, is nontrivial.

To answer fundamental questions about spatiotemporal compartmentalization of
cellular signaling, building on advances in meshing, we introduce Spatial
Modeling Algorithms for Reaction and Transport (SMART), a Python-based
software package to construct and solve systems of mixed-dimensional reaction-
transport equations [24, 25]. SMART complements other software packages for
computational cell biology such as VCell [26, 27] and MCell [28] by providing
a unique toolkit to solve spatial signaling networks over complicated
geometries using the finite element method (FEM). Similar to systems biology
frameworks such as the Systems Biology Markup Language (SBML) [29], SMART
takes user input that specifies species, reactions, compartments, and
parameters involved in a given biological network. Herein, we describe this
workflow and then demonstrate the use of SMART in several biological test
cases, including YAP/TAZ-mediated mechanotransduction in cells on
micropatterned substrates, calcium signaling in neuronal dendritic spines,
calcium release within a cardiomyocyte, and synthesis of ATP in a
mitochondrion. Several of these examples use realistic subcellular geometries
derived from electron micrographs and conditioned with GAMer2. We further
demonstrate the accuracy of SMART via a series of addition numerical
verification tests and summarize its performance and scalability.

## Results

Figure 1:  Mixed dimensional reaction-transport networks in cells. A)
Schematic of reaction-transport system in SMART. Given the topological
relationships between different compartments in cells and information on
reactions between species, fluxes across boundaries, and diffusion rates of
individual species, SMART assembles a finite element system of equations.
Panel created with Biorender.com. A single model in SMART may include both
volume species (circles in inset) and surface species (rectangles in inset)
that can all diffuse and react with one another. B) Whole cell geometry with
segmented organelles from volume electron microscopy. Adapted from Figure 1a
in Heinrich et al. 2021 [8], permission pending. C) Calcium release from the
endoplasmic reticulum (ER) in a realistic geometry. For illustration, a linear
molecular leak flux from the ER ($j_{leak}=j_{0}(c_{ER}-c_{i})$) was assumed
starting at $t=0$. This realistic ER geometry was derived from electron
micrographs of dendritic spines [30].

### Spatial modeling algorithms for reaction and transport to simulate
biological networks

Cell signaling networks rely on non-trivial chains of molecular reactions and
transport mechanisms acting within and across compartments: extracellular,
intracellular, and subcellular spaces, and membranes (Figure 1A). These
cellular domains define three-dimensional bulk _volumes_ while the membranes
can be viewed as two-dimensional manifold _surfaces_. The geometry of these
compartments can be accurately represented through synthetic or imaging-guided
computational meshes in the form of simplicial complexes, with individual
volumes and surfaces identified and labeled by tags [19] (Figure 1B–C,
Methods). Spatial modeling of cellular pathways and processes describe the
distribution and evolution of different _species_ present in or on these
compartments; e.g., ion concentrations such as Na^+ in the cytosol or
interstitial fluid, Ca^2+ in subcellular compartments, or the density of
receptors or other channels distributed along the plasma membranes.

Mathematically, we describe diffusion of such species and reactions between
species, within or across compartments, via coupled systems of time-dependent,
nonlinear and mixed-dimensional partial differential equations defined over
the computational geometries (Figure 2, Methods). Our modeling assumptions
enable species to diffuse within volumes and on surfaces, and to be
transported across surfaces to cross between volumes. Reactions between a
single or multiple species occur within compartments (volume or surface
reactions) or in adjacent compartments (volume-surface or volume-surface-
volume reactions) (Figure 2A). The use of a mixed finite element
discretization in space allows for high numerical accuracy and geometric
flexibility [16]. Crucially, this approach allows for spatial variation of
each species within compartments and transport across compartments, while
conserving mass and momentum. The computational model is allowed to evolve
over time, yielding detailed predictions for changes in each species within
the specified geometries (Figure 2E). High performance sparse numerical linear
algebra [31] in combination with scalable finite element algorithms [23]
allows for computational models with millions of degrees of freedom to be
solved efficiently.

In the following sections, we demonstrate SMART via a series of use cases that
model different signaling networks in cellular and subcellular geometries. Our
examples begin at the whole-cell scale and approach the single organelle
scale, including both idealized geometries as well as experimentally derived
geometries from 3D electron microscopy.

Figure 2:  SMART workflow. A) Illustration of the basic components used to
define a model in SMART - species, reactions, parameters, and compartments.
The graphic contains 3 volume compartments ($\Omega^{m}$) and 3 surface
compartments ($\Gamma^{q}$), with illustrations of the reaction types
supported by SMART, including volume, surface, volume-surface, and volume-
surface-volume. Panel created with Biorender.com. B) Model geometry specified
by mesh generated in Gmsh or conditioned in GAMer2. C) Assembly of equations
from reaction specifications in SMART. Each volume species and each surface
species has an associated PDE with boundary conditions, as shown on the right
(equations in Section A.3). D) SMART model discretization using finite
elements. The nonlinear system is discretized using linear finite elements in
FEniCS and then the block matrix problem is assembled in PETSc. The matrix is
nested in terms of compartments involved in each interaction, as summarized on
the right. E) Model solution and postprocessing in SMART. The system is solved
iteratively at each time step until reaching $t_{final}$. Results are post-
processed to examine changes in concentration over time and space.

### YAP/TAZ mechanotransduction in cells on micropatterned substrates

To demonstrate use of SMART at the whole-cell length scale, we consider the
model of YAP/TAZ mechanotransduction originally developed by Sun et al. [32]
and extended to a spatial model by Scott et al. [33]. This model considers the
intracellular signaling cascade induced by a cell adhering to a substrate with
a given stiffness, from phosphorylation of focal adhesion kinase (FAK) to
downstream activation of the actomyosin cytoskeleton and nuclear translocation
of the transcriptional regulatory proteins Yes-associated protein (YAP) and
PDZ-binding motif (TAZ) (Figure 3A). Generally, on stiffer substrates,
increased FAK phosphorylation results in increased actin polymerization and
myosin contractility, leading to dephosphorylation of cytosolic YAP/TAZ and
subsequent translocation of YAP/TAZ into the nucleus.

Here, we test the predicted effects of cell adhesion to micropatterned
surfaces on the localization of YAP/TAZ to the nucleus. We consider three
different contact region geometries previously tested experimentally [34] \-
circular, rectangular, and star-shaped micropatterns (Figure 3B-C). In all
cases, the cell volume and size of the contact region are the same, but the
total plasma membrane surface area is increased for geometries with greater
curvature at the contact regions (see Appendix E for details on generating
these meshes).

Figure 3: Spatial model of YAP/TAZ mechanotransduction for cells on
micropatterned surfaces. (cont. on next page) Figure 3: (cont. from previous
page.) A) Schematic of YAP/TAZ mechanotransduction signaling pathway.
Phosphorylation of FAK within the region of cell-substrate adhesion leads to
cytoskeletal activation, triggering the dephosphorylation of YAP/TAZ and the
opening of nuclear pore complexes (inset), allowing for the transport of
YAP/TAZ into the nucleus. Panel created with Biorender.com. B) Measurements of
cytoskeletal activation in cells on micropatterned substrates. Modified from
Figures 2 and 4 in [34], permission pending. C) Summary of geometries for
cells spread on circular, rectangular, and star-shaped micropatterns. D)
Simulations of YAP/TAZ mechanotransduction in cells on circular, rectangular,
and star-shaped micropatterns. For each case, the side view is pictured,
including both F-actin in the cytosol and YAP/TAZ in the nucleus, as well as
the bottom-up contact region view showing local activation of actin
polymerization in regions of higher curvature along the cell contour. Bottom-
up views also include contours showing lines of constant concentration. E-F)
F-actin (E) and nuclear YAP/TAZ (F) dynamics plotted over time for all three
cases.

We consider the predictions of this model across all three geometries on a
very stiff substrate such as glass. We find that regions of the cell where the
plasma membrane surface area to cytosolic volume ratio is the highest have
high concentrations of signaling molecules of interest. In particular, we
observe elevated levels of F-actin in these regions, in agreement with the
results from Maduram et al. [34] (Figure 3B, D, Videos 1, 2 and 3).
Micropatterns with highly curved regions along the perimeter (star or
rectangular patterns) show greater increases in overall cytoskeletal
activation (Figure 3E) and, consequently, higher nuclear abundance of YAP/TAZ
(Figure 3F). Importantly, regardless of the shape of the contact area, all
increases in YAP/TAZ are much lower than those predicted by a well-mixed model
of YAP/TAZ mechanotransduction (compare to Figure 6D) due to persistent
gradients in cytoskeletal activation over the cellular geometry. That is, the
F-actin concentration at the nuclear membrane is much lower in this spatial
model compared to the well-mixed case, in which the effects of signal
attenuation in regions of the cytosol further away from the plasma membrane
are neglected. This observation highlights the importance of accounting for
spatial effects in cell signaling networks.

### Intracellular calcium dynamics in realistic subcellular geometries

We next consider spatial models defined over realistic _subcellular_
geometries acquired from 3D electron microscopy. We utilize previously
published models of calcium dynamics within dendritic spines in neurons [35]
and cardiomyocyte calcium release units (CRUs) [12]. Each case includes
calcium influx through the plasma membrane and calcium exchange across the
sarcoplasmic or endoplasmic reticulum (SR or ER) membrane, as well as calcium
binding and unbinding to buffering proteins within the cytosol and SR/ER
(Figure 4A). The calcium dynamics in each system are influenced by the
geometry and relative spatial arrangement of organelles [36, 37].

Using the model previously implemented for idealized dendritic spine
geometries by Bell et al. [35], we simulate calcium changes within a realistic
dendritic spine containing a specialized form of endoplasmic reticulum known
as the spine apparatus (SA) [30]. Calcium influx occurs through voltage-
sensitive calcium channels (VSCCs) located in the head and a section of the
neck and through NMDA receptors localized to a dense region of proteins known
as the post-synaptic density (PSD) (Figure 4B). Each of these fluxes is
written as an analytical expression over time, and dependent on a specified
change in the membrane potential. Calcium is removed from the cytosol via
efflux across the plasma membrane through the sodium calcium exchanger (NCX)
and plasma membrane calcium ATPase (PMCA) and influx into the spine apparatus
through the sarco/endoplasmic calcium ATPase (SERCA). Over the short time
scale considered in this model, we assume that calcium entry into the SA
dominates over release. Simulations reveal that calcium elevation is highest
within the spine head, also leading to a more substantial increase in calcium
concentration in the SA located in this region (Figure 4C-D, Video 4). The
peak calcium approaches 5 $\upmu$M, lower than the 8 $\upmu$M peak observed in
the original model [35]. This lower peak calcium concentration is accounted
for by a higher surface area of SA present in the spine head for this
realistic geometry compared to the smaller SA in the original simulations
(here: 10.6 $\upmu$m2 SA per $\upmu$m3 cytosol; original: 2.31 um2 SA per um3
cytosol).

Given the importance of calcium in signaling pathways, to ensure robustness
across different model formulations and geometries, we also use SMART to model
calcium release from the SR in a CRU. The CRU geometry[38] includes a
continuous section of SR, one T-tubule volume, and two mitochondria. We do not
consider any species within the mitochondria or T-tubule, but rather assume
that the T-tubule calcium concentration remains constant due to its continuity
with the extracellular space and treat the mitochondria as passive diffusion
barriers, in line with previous modelling by Hake et al.[12]. Similarly, the
NCX, PMCA, and leak fluxes through the T-tubule membrane, the ryanodine
receptor (RyR) and SERCA fluxes through the SR membrane, as well as calcium
buffering due to calmodulin, troponin, and ATP in the cytosol and
calsequestrin in the SR, all follow established relationships[12]. However, we
here utilize a full spatial discretization of the sarcoplasmic reticulum
interior, which was previously treated as a series of well-mixed
subregions[12]. For comparison, we consider two conditions – one with SERCA
present throughout the SR membrane and another with no SERCA activity. In
agreement with the original model, calcium reaches a concentration of several
$\upmu$M near the T-tubule / SR junction and about 1 $\upmu$M further away
from the sites of RyR release (Figure 4E-F, Video 5). This calcium spike is
short-lived, as the model assumes that RyRs close upon sufficient reduction of
SR calcium. Including active SERCA in the SR membrane results in a slightly
prolonged calcium elevation and robust refilling of the SR calcium store over
time (Figure 4F, Video 6), again in line with previous findings[12]. In both
calcium signaling examples presented here, the spatial behavior captured by
SMART is critical to the overall dynamics, with the proximity between calcium
sources and sinks determining the overall extent of calcium increase.

Figure 4:  Calcium dynamics within a realistic dendritic spine and calcium
release unit. (cont. on next page)  Figure 4: (cont. from previous page.) A)
General schematic of a calcium signaling cascade in a subcellular region. In
response to biochemical or electrical stimulus, calcium influx occurs through
the PM and calcium is released from the ER/SR store. Calcium elevations are
counteracted by calcium efflux into the extracellular space and repackaging
into the ER/SR through SERCA. Finally, the changes in free calcium are
dampened due to calcium binding to proteins in the cytosol (buffering). Panel
created with Biorender.com. B) Realistic geometries of a dendritic spine and a
cardiomyocyte calcium release unit derived from electron microscopy. The
dendritic spine branches off from the main dendritic shaft and contains a
lamellar ER structure known as the spine apparatus (SA, gold), as well as a
denser region of signaling proteins in a subregion of the spine head, the
postsynaptic density (PSD, red). The calcium release unit consists of a
section of SR (gold) closely apposed to T-tubules (dark grey) in a region
known as the junctional cleft. Mitochondria (light grey) serve as passive
diffusion barriers in this case. C) Simulation of calcium dynamics in a
dendritic spine. Calcium elevations are concentrated to the dendritic head;
following the initial influx, calcium is pumped into the SA via SERCA. D)
Plots of cytosolic and SA calcium dynamics in the dendritic spine head vs. the
dendritic shaft. E) Simulation of calcium dynamics in a CRU. Calcium release
from the SR results in a calcium spike near the center of the geometry,
followed by refilling of the SR with calcium. F) Plots comparing the dynamics
of cytosolic and SR calcium with vs. without SERCA.

### ATP synthesis in realistic mitochondrial geometries

We finally consider a biochemical network within a single organelle, modeling
the generation and transport of ATP within a realistic mitochondrial geometry.
This mitochondrial geometry was previously reconstructed from serial electron
tomograms and the resulting mesh was conditioned in GAMer2 [14, 39]. It was
then used for simulations of ATP generation in a well-mixed context and for
particle-based simulations in MCell [14, 40]. Here, we implement a
thermodynamically consistent continuum-based model version[40] and compare our
current results to those from well-mixed ODE simulations and particle-based
simulations[14, 40].

Figure 5:  ATP synthesis in a realistic mitochondrial geometry. A) Schematic
of the ATP synthesis process in mitochondrion. ADP in the mitochondrial matrix
is converted into ATP by ATP synthase in the cristae membrane, then
transported into the inner membrane space through ANTs and exported into the
cytosol through VDACs. Panel created with Biorender.com. B) Realistic
mitchondrial geometry. The inner membrane divides two volume compartments -
matrix (yellow) and inner membrane space (blue). The inner membrane itself is
divided into a portion on the periphery (yellow) and invaginations known as
cristae (brown). C) Dynamics of IMS ATP and matrix ADP in a mitochondrion with
uniform distributions of ATP synthase and ANTs in the inner membrane. D-E) IMS
ATP concentration (D) and cytosolic ATP concentration (E) over time for
uniform vs. cristae-localized distributions of IMS proteins. ODE model
predictions are plotted as well for ease of comparison (dashed lines).

The model considers the generation of ATP in the mitochondrial matrix by ATP
synthase, followed by the transport of ATP into the inner membrane space (IMS)
by adenine nucleotide transporters (ANTs) and export into the cytosol through
voltage-dependent anion channels (VDACs) in the outer membrane (OM) (Figure
5A-B). Previous experimental data and modeling results suggest that the
spatial organization of ATP synthase and ANTs is a key determinant of the
magnitude of ATP changes in the cytosol [14]. Accordingly, we tested two
alternative spatial arrangements of ATP synthase and ANTs while keeping the
total number of molecules unchanged. In the first case, we distributed these
molecules uniformly throughout the inner membrane (IM), whereas in the second,
more physiological case, we colocalized ANTs and ATP synthase in invaginations
of the IM known as cristae. In both cases, we maintained the same total amount
of ANTs and ATP synthase molecules.

A uniform distribution of ANTs results in dynamics of IMS and cytosolic ATP
similar to those predicted by the ODE model and particle-based
simulations(Figure 5D-E). ATP in the IMS initially decreases rapidly as it
binds to ANTs, then gradually increases as ATP synthase produces more ATP. In
comparison, concentrating ANTs and ATP synthase in the cristae results in a
less pronounced initial reduction of ATP in the cytosol. This phenomenon was
previously termed “energy buffering” [14] and was attributed to the spatial
separation between ANTs in the cristae and VDACs in the OM. Rapid changes to
ATP concentration in the inner cristae space do not immediately affect ATP
levels in the cytosol as time is required for diffusion between the outer
membrane and cristae. This delay and dampening effect could make the cell more
resilient to a noisy environment featuring rapid changes in the availability
of ATP and ADP. This spatial phenomenon is captured via SMART simulations but
is not accessible via well-mixed ODE models.

### Verification and validation of the SMART computational algorithms.

To verify our numerical and computational approach, we examine the simulation
results for an example with an analytical solution at steady state (details in
Section D.1). In this example, originally examined by Meyers et al. [41], a
protein is phosphorylated at the cell membrane and dephosphorylated throughout
the cytosol (Figure 6A). We considered the case of a thin slab whose spatial
solution is well-described by the 1D solution along its thickness. This
geometry mimics the case of an adherent cell spread out on a flat substrate.
To demonstrate spatial convergence, we uniformly refined the reference mesh
four times, resulting in maximum cell sizes of 3 $\upmu$m, 1.5 $\upmu$m, 0.75
$\upmu$m, and 0.375 $\upmu$m. Similarly, to show temporal convergence we
halved the time step 7 times from 0.64 s to 0.01 s. Figure 6B-C show the
temporal and spatial convergence respectively, for three different diffusion
coefficients, using the finest mesh (B) and finest time step (C). We observe
that the numerical error is consistently lower for higher diffusion
coefficients. Moreover, the mesh resolution error dominates up to some
critical value of the time step for each $D$, as indicated by the plateau in
each curve (Figure 6B). This critical time-step is smaller for high values of
$D$, reflecting the need for smaller time-steps to achieve optimal convergence
in cases of rapid diffusion. At the smallest time step, mesh refinement error
dominates in all cases, as demonstrated by the theoretically-expected and
optimal second-order convergence in the $\mathcal{L}_{2}$ norm of the error
with respect to element size $h$ (Figure 6C).

Next, we focus on the accuracy and convergence of the numerical results for
our three main biological scenarios. In the case of YAP/TAZ
mechanotransduction, we can compare our numerical solutions to solutions of
the well-mixed ODE approximation. Considering the case of a cell with a
circular contact region on a glass substrate, we set diffusion coefficients
for all volume species to 1000 and treat all surface species as uniform and
non-diffusing. We start with a time step of $\Delta t=0.01$ s and adaptively
adjust the step according to the number of Newton iterations required for
convergence at the previous time (see Table S1 for details of our adaptive
time-stepping scheme). Upon mesh refinement, the average F-actin concentration
and nuclear YAP/TAZ concentration converge to a solution close to the
trajectories predicted by the ODEs (Figure 6D). For the main downstream
quantity of interest, YAP/TAZ nuclear concentration, the refined mesh differs
a maximum of 0.6% from the well-mixed solution. Note that the most refined
mesh here is coarser than the mesh used in simulations shown in Figure 3 (see
all mesh statistics in Table S9 and LABEL:tab:mech-compartments).

Dendritic spine simulations were initially conducted on a mesh with
$h_{min}=0.009$ $\upmu$m and $h_{max}=0.199$ $\upmu$m and a time step of
$\Delta t=0.00025$ s. After two refinements of the mesh, we observe that the
average calcium concentrations are similar across refinements, but that the
solutions exhibit local spatial differences (Figure 6E). In general, lower
peak values are observed in the more refined simulation. However, across
refinements, the peak calcium in the spine head deviated by a maximum of
0.52%, demonstrating the accuracy of the baseline results (Figure 6E).
Similarly, testing time steps from 0.00025 s to 0.002 s reveals convergence to
a single solution consistent with that in Figure 4 (Figure S2).

Finally, we assess the accuracy of our solutions at the single organelle level
by comparing our mitochondrial ATP production simulations to a well-mixed
approximation. The mitochondrion mesh resolution is $h_{min}=0.00661$ $\upmu$m
and $h_{max}=0.0383$ $\upmu$m, and our time step is adaptively adjusted from
an initial value of $\Delta t=0.02$ ms. When setting the diffusion
coefficients for all volumetric nucleotide species to 150, we obtain a close
match to the ODE solution (Figure 6F). In particular, the predicted value of
cytosolic ATP deviates by a maximum of 0.2% from the well-mixed solution
(Figure 6F).

Figure 6: Numerical verification and validation of SMART. A) Model of protein
phosphorylation at the plasma membrane. B-C) Convergence of system describing
phosphorylation of proteins at the plasma membrane (as described in Ref. 41)
upon time step refinement (B) and mesh refinement (C). The $\mathcal{L}_{2}$
error between the computed and analytical solutions for all combinations of
temporal and spatial refinement are shown in Table S6, Table S7 and Table S8
for diffusion coefficients $D$ 10, 100, and 1000 respectively. D) Comparison
of ODE solution to full SMART simulation for mechanotransduction example with
fast diffusion at different mesh refinements. E) Spatial differences in
solution for dendritic spine calcium for course vs. refined mesh. Maximum and
average calcium concentrations in different regions are reported for the
coarse vs. fine mesh. F) Comparison of ODE solution to full SMART simulation
of mitochondrial ATP production for fast-diffusing nucleotide species.

### Computational performance and scalability

To study the scalability of our computational framework, we consider the
simulation of intracellular calcium dynamics in an electron micrograph-based
representation of a dendritic spine (Example 2). Starting with the baseline
computational mesh, we consider two additional levels of uniform mesh
refinement resulting in three discrete representations (‘standard’, ‘fine’,
and ‘extra fine’). These meshes are composed of $106\,979$ ($18\,649$),
$855\,832$ ($146\,024$), and $6\,846\,656$ ($1\,154\,871$) tetrahedral cells
(vertices), respectively, thus corresponding to an $8\times$ increase in the
number of cells for each refinement and a similar increase ($7.8-7.9\times$)
in the number of vertices. The corresponding numbers of degrees of freedom
($N$, dimensions of the discrete solutions) are: $49\,194$, $323\,328$,
$2\,299\,090$; thus corresponding to a $6.6\times$ and $7.1\times$ increase
between refinements. For each simulation, we use a 0.001 s timestep through
$t=0.025$ s. We first note that the number of nonlinear iterations per time
step is constant ($4$) across time steps and refinement levels. Next,
inspecting the total run time of each simulation as a function of the
computational cost (measured in terms of the number of degrees of freedom
$N$), we observe that the simulation time scales close to log-linearly
(between $\mathcal{O}(N)$ and $\mathcal{O}(N\log N)$) (Figure 7A). The
computational cost associated with the initialization of the SMART symbolic
problem representation is small ($2.0\%$, $2.4\%$, $3.6\%$). A break-down of
the total run times shows that the simulation time is persistently dominated
by the finite element assembly ($87.4\%$, $88.4\%$, $87.0\%$), followed by the
iterative solution of the linear systems ($9.5\%$, $8.4\%$, $8.8\%$) (Figure
7B).

Figure 7: Performance and scalability of SMART for dendritic spine calcium
dynamics A) Total run times for 25 simulated time steps of intracellular
calcium dynamics for increasing mesh refinement levels (standard, fine and
extra fine), averaged over 5 runs for standard and fine meshes and 3 runs for
the extra fine mesh. B) Relative contributions to the total run time for each
level, averaged over 5 runs for standard and fine meshes and 3 runs for the
extra fine mesh. The computational cost is dominated by the finite element
assembly and function evaluation due to the complexity of the coupled
nonlinear equations.

## Discussion

Cell function is tightly linked to cell shape, as evidenced by the diverse
shapes exhibited by different cell types, from neurons with branch-like
extensions to cardiomyocytes that assemble in block-like morphologies [42, 43,
44]. The importance of geometry in function extends to the single organelle-
level; for instance, the inner membrane of mitochondria forms tortuous
invaginations (cristae) that facilitate efficient production of ATP. The
importance of cell and organelle geometry are generally well-appreciated, but
detailed geometries have proven difficult to include in models of cell
signaling. The models showcased within this work exhibit the use of our
simulation technology and software, SMART, to simulate such signaling networks
with spatial resolution, with a particular focus on realistic cell and
organelle geometries.

SMART offers many possibilities for exploring the impacts of cell shape in
spatiotemporal signaling models in cell biology. This opportunity is driven by
a wealth of new imaging data from modalities such as volume electron
microscopy and super-resolution fluorescence microscopy [8, 9, 10, 11]. While
the experimentally-derived cell geometries here were all extracted from
electron micrographs, other types such as e.g. fluorescence data could be
utilized for future applications of SMART. Fluorescence data could offer the
additional opportunity to inform not just cell and organelle geometry, but
also the spatial localization and dynamics of different molecules such as
membrane receptors. There are a variety of high-quality public datasets
available from super-resolution microscopy and volume electron microscopy
alike due to projects such as the Allen Institute’s Cell Explorer [45, 46] and
Janelia’s OpenOrganelle [8, 47]. Biophysical modelling leveraging such imaging
datasets takes full advantage of the recent efforts to improve mesh
conditioning for biological samples such as GAMer2 [19], VolRover [48] or
fTetWild [49].

Indeed, the examples in this paper demonstrate the fundamental importance of
geometry in models of cell signaling. For instance, in our model of YAP/TAZ
mechanotransduction, while results agree qualitatively between a well-mixed
model and a full spatial description, quantitative predictions of the model
differ considerably between these two modeling regimes. Furthermore, certain
aspects of the signaling network cannot be captured by a well-mixed model,
such as increased levels of actin polymerization near the plasma membrane
compared to elsewhere in the cell. Similarly, despite the fast diffusion of
small molecular species such as calcium and ATP, we find that spatial factors
such as the adjacency between sources and sinks (dendritic spine and CRU
examples) or the distribution of surface species (ATP example) are key
determinants of system behavior. Indeed, in the case of ATP generation,
changing the distribution of ATP synthase and ANTs in the inner membrane has a
non-negligible effect on ATP dynamics in the cytosol. Considering that model
parameters are commonly estimated using well-mixed models, parameterizing
spatial models remains an important challenge. SMART is well-positioned to
address these issues via auxiliary tools such as dolfin-adjoint for
sensitivity analysis and parameter estimation in FEniCS-based models [50].

In addition to SMART, there exist several complementary software options
available to computational biologists and biophysicists for examining the
spatiotemporal behavior of cell signaling networks. Among the most popular are
Virtual Cell (VCell) [26, 27] and MCell [28], both open source projects
focused on continuum and particle-based models of biochemical networks in
cells, respectively. Within this broader context, SMART offers unique
capabilities and features for those users wishing to test complex cell
geometries in a continuum framework. For instance, VCell provides a robust
option for solving mixed dimensional reaction-diffusion equations using the
finite volume method, but currently only pixelized or voxelized grid meshes
are supported [51]. MCell, in contrast, does support surfaces defined by
unstructured meshes; however, MCell uses particle-based simulations to
describe stochastic reaction-diffusion networks. Such simulations are much
more accurate when considering a small number of molecules but cannot yet be
scaled to whole-cell simulations where many billions of molecules may be
present. Alternative stochastic simulation packages such as the STochastic
Engine for Pathway Simulation (STEPS) [52, 53] or Lattice Microbes [54] use
other algorithms to model stochastic reaction-diffusion networks within cells.
Many other options allow users to readily assemble signaling networks and/or
perform parameter estimation and sensitivity analysis in the well-mixed case
[55, 56]. These complementary approaches work well alongside SMART to help
develop robust models of cell signaling and select parameters for spatial and
well-mixed models alike.

In terms of limitations, SMART currently only supports describing a set of
compartments of codimension 1, i.e., only 3D/2D or 2D/1D, whereas in
principle, 1D or even 0D features could contribute in the full 3D case. We
note that 0D/well-mixed features can currently be included by explicitly
updating parameters at each time step, as in our ATP generation example (see
Equations 1 and 2). SMART also currently only supports diffusion as a
transport mechanism, whereas our framework can easily be expanded in the
future to support e.g. reaction-advection-diffusion problems. This can also
naturally be extended to the case in which the cell boundary moves over time
according to some prescribed velocity field. SMART does not currently consider
the effects of other physical processes, such as electrodiffusion or
mechanical forces, that are often tightly coupled to signaling. In fact, this
remains an issue across most of computational cell biophysics. FEM-based
approaches, such as that of SMART, offer one promising approach to efficiently
solve such problems, due to their general applicability to solve a range of
PDEs. In summary, SMART represents an advancement towards realistically
simulating cellular and subcellular systems. The ability to do so will not
only drive forward the understanding of cell biology, but could act as a tool
to predict the response of cells to diverse treatments at an unprecedented
level of detail.

## Methods

The mathematical models, computational algorithms and simulation technology
underlying SMART are summarized here. An extended description is provided in
Appendix A.

### Domains, geometries and interfaces

We consider model geometries embedded in 2D or 3D represented by tessellated
curves, surfaces and volumes. Each geometry is described by a collection of
three-dimensional volumes (or two-dimensional surfaces in the 2D case) with
boundaries and interfaces between volumes represented by two-dimensional
surfaces (or one-dimensional curves in the 2D case). We refer to the highest
dimensional compartments as the _bulk domains_ ($\Omega$) and the lower
dimensional boundaries or interfaces as the _surfaces_ ($\Gamma$). Moreover,
each of these domains are represented by a simplicial tessellation formed by
the mesh cells (intervals, triangles or tetrahedra), mesh facets (points,
intervals or triangles), and mesh vertices (points). Importantly, all
(sub)regions, boundaries and interfaces are defined relative to a single
overarching _parent_ mesh [16]. Subdomains $\Omega^{m}\subseteq\Omega$ are
defined as a set of mesh cells with a common label or tag $m\in\mathcal{M}$,
and similarly exterior boundaries and interior interfaces $\Gamma^{q}$ are
defined as a set of mesh facets again with a common tag $q\in\mathcal{Q}$. In
the examples presented here, we generated parent meshes and labeled their
subdomains using either GAMer2 [19] to construct meshes from electron
microscopy data or Gmsh [57] to generate idealized cell geometries.

### Coupled multi-domain reaction-transport equations

SMART represents the spatially and temporally-varying concentrations of
multiple species coexisting within domains via a general abstract modelling
framework based on first principles. Each concentration is defined over a
subdomain $\Omega^{m}$ and/or surface $\Gamma^{q}$. The evolution and
distribution of these concentrations are described by a coupled system of
time-dependent and nonlinear partial differential equations representing
conservation of mass, the diffusion of each species, reactions between species
within the bulk or surface domains, and fluxes between the bulk and surface
domains (see Section A.3 for a full description of the general framework).
Reactions are assumed to be local in the sense that different species may
interact within each subdomain $\Omega^{m}$, and on or across the surfaces
$\Gamma^{q}$ via the surface itself and its neighboring subdomains. Specific
computational models specify the effective diffusivity of each species and
symbolic expressions for model reactions and fluxes, the initial concentration
of each species, and any bulk or surface source or sinks. Input parameters may
be constant or spatially or temporally varying. For further details, we refer
to Section A.4.

### Numerical approximation and solution strategies

Each system of mixed-dimensional reaction-transport equations is discretized
in time using a first-order accurate implicit Euler scheme with a uniform,
variable or adaptive timestep size, yielding a coupled system of nonlinear
partial differential equations at each time step. For the spatial
discretization, we employ a monolithic finite element method via the FEniCS
finite element software suite [23]. As unknown discrete fields, we consider
the concentration of each species $u_{i}^{m}$ defined in each bulk subdomain
$\Omega^{m}$ and the concentration of each species $v_{j}^{q}$ defined over
each surface $\Gamma^{q}$. All discrete fields are represented by continuous
piecewise linear finite elements defined relative to the mesh used to
represent the geometry. Since the equations are coupled across fields and
domains, the resulting nonlinear systems of discrete equations take a block
structure [16] (Figure 2D).

All nonlinear systems of equations are by default solved by Newton-Raphson
iterations with an exact, automatically- and symbolically-derived Jacobian.
The time step for each system is either set as uniform throughout the course
of the simulation or can be set adaptively based on the number of Newton-
Raphson iterations required for convergence at the previous time step, as
summarized in Table S1. In case of nonlinear solver divergence or negative
solutions, repeated restarts with associated reductions in time step are
invoked as an optional mediation strategy. The linear systems are solved
iteratively using Krylov solvers in PETSc [31]. For default solver tolerances
and settings, we refer to the open-source SMART code [25].

### SMART model specifications

As outlined in the first results section, the user must specify species,
reactions, parameters, and compartments, and link these to a parent mesh prior
to initializing a simulation. Several minimal use cases are given in the SMART
documentation [25] and the code for all examples in this paper is freely
available [58]. In short, each species, reaction, parameter, or compartment is
defined as a Python object, each with associated properties detailed in
Appendix C. All instances of a given object type are then stored in an ordered
Python dictionary, resulting in a single “container” for each type. The parent
mesh is either generated using Gmsh or read directly from an hdf5 or xml file.
The containers and mesh are then used to initialize the SMART model.

### Model setup for biological test cases

Here, we briefly outline some of the relevant details for each test case; for
full model specifications, refer to Appendix F and our code [58].

#### Example 1

This example is described in a previous publication that conducted similar
simulations in VCell [33]. There are 24 species in the original model, 11 of
which are eliminated after accounting for mass conservation. We note that this
simplification is only possible when different forms of a given molecule
reside in the same compartment and share the same diffusion coefficient. This
simplification, for instance, cannot be applied to actin species (F-actin and
G-actin), as they have drastically different diffusion coefficients. All
parameters were used unaltered from Scott et al.[33], with the exception of
the altered diffusion coefficients for well-mixed simulations.

We generate meshes for these cell geometries as described in Appendix E. We
minimized computational expense by exploiting the symmetries of each geometry;
in the case of a circular contact region, the geometry is axially symmetric,
allowing us to simulate the model over a 2-dimensional mesh while ensuring the
correct $r$-dependence in all integrals (Section C.3). The rectangular contact
region has two symmetry axes, allowing us to only simulate one-quarter of the
full geometry, treating the faces on the symmetry axes as no flux boundaries.
Similarly, the star contact region has five symmetry axes, allowing us to
simulate one-tenth of the full mesh. Simulations were run to $t=10,000$ s, by
which time the system has plateaued to an apparent steady-state.

#### Example 2

The model of calcium dynamics within a dendritic spine was derived from Bell
et al.[35]. This model involves only four dynamical species - calcium, fixed
calcium buffers bound to the plasma membrane, mobile calcium buffers in the
cytosol, and calcium in the spine apparatus, which was previously assumed to
be constant. Calcium bound to buffers was not treated as a separate species,
but was implicitly included assuming mass conservation and the same diffusion
coefficient for buffering protein and buffering protein bound to calcium. All
other quantities that change over time were treated as time-dependent
parameters, explicitly defined by functions of time. As in the original model,
NMDA receptors were restricted to the postsynaptic density. VSCCs were
restricted to the spine head and a portion of the neck to match the region of
stimulus in Bell et al. [35]. All other membrane fluxes were uniform
throughout the spine head, neck, and dendritic shaft.

Our analysis of the calcium release unit used the model presented by Hake et
al. [12] with minor modifications. In their simulations, the SR was treated as
a composite of several well-mixed compartments and the wider space surrounding
the CRU was simulated as connected well-mixed compartments. Instead, we
included the entire SR volume explictly in our spatial simulations, and
boundaries of the cytosolic mesh were treated as no-flux surfaces. As above,
complexes between calcium and buffering species (ATP, CMDN, TRPN, and CSQN)
were implicitly considered under the assumptions of buffer mass conservation
and unchanged diffusion coefficients. The main mechanism for calcium
elevations in this model is release from the SR through ryanodine receptors,
which is assumed to terminate when calcium concentration falls below a certain
threshold. Rather than explicitly encoding this discontinuity in the
equations, we enforce this condition manually–we assign zero conductance to
the RyRs after the average calcium concentration in the entire SR reaches the
threshold. SERCA channels were either included uniformly through the SR
membrane or were excluded entirely.

#### Example 3

The model of ATP generation was directly adapted from the model by Garcia et
al. [40]. The model considers 6 states of ATP synthase and 9 states of ANTs,
as well as ATP concentration in the matrix and inner membrane space (IMS) and
ADP concentration in the matrix. Assuming mass conservation, only 5 states and
8 states need to be modeled explicitly for ATP synthase and ANTs. ADP
concentration in the IMS is assumed to be constant and ATP concentration in
the cytosol was solved for using the following coupling scheme. Prior to
solving the PDEs at each time step ($t=t_{n}$), the ATP concentration
$T_{cyto}$ at the next time point $t_{n}+\tau_{n}$ was estimated as:

$T_{cyto,est}(t_{n}+\tau_{n})=T_{cyto}(t_{n})+\tau_{n}\frac{10^{18}\text{mM}\upmu\text{m}^{3}/\text{mol}}{vol_{cyto}N_{A}}\iint\limits_{\Gamma_{OM}}k_{vdac}[VDAC](T_{IMS}(t_{n})-T_{cyto}(t_{n}))d\Gamma,$
(1)

where $T_{IMS}$ is the ATP concentration in the IMS, $vol_{cyto}$ is the
volume of cytosol immediately surrounding the mitochondrion, $N_{A}$ is
Avogadro’s number, and other parameters are summarized in Appendix F. This
updated value of $T_{cyto}$ was then used in solving the PDE system, after
which the estimated value was updated for consistency with the implicit Euler
time discretization:

$T_{cyto}(t_{n}+\tau_{n})=T_{cyto}(t_{n})+\tau_{n}\frac{10^{18}\text{mM}\upmu\text{m}^{3}/\text{mol}}{vol_{cyto}N_{A}}\iint\limits_{\Gamma_{OM}}k_{vdac}[VDAC](T_{IMS}(t_{n}+\tau_{n})-T_{cyto,est}(t_{n}+\tau_{n}))d\Gamma.$
(2)

## Acknowledgements

We would like to acknowledge the members of the Rangamani Lab for their
valuable insights throughout the development of SMART. We also thank Dr.
Cécile Daversin-Catty for her guidance on using mixed dimensional features in
FEniCS.

E.A.F. was supported by the National Science Foundation under Grant #
EEC-2127509 to the American Society for Engineering Education and funding from
the Wu Tsai Human Performance Alliance at UCSD to P.R. P.R. was supported by
the Air Force Office of Scientific Research (AFOSR,
https://www.wpafb.af.mil/afrl/afosr/) Multidisciplinary University Research
Initiative (MURI) FA9550-18-1-0051. J.G.L. was supported by a fellowship from
the UCSD Center for Transscale Structural Biology and Biophysics/Virtual
Molecular Cell Consortium and C.T.L. was supported by a Hartwell Foundation
Postdoctoral Fellowship and a Kavli Institute of Brain and Mind Postdoctoral
Fellowship.

M.E.R was supported by the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme under grant agreement
714892 (Waterscales), by the Research Council of Norway under grant #324239
(EMIx), and by the Foundation Kristian Gerhard Jebsen via the K. G. Jebsen
Center for Brain Fluid Research, and by the Fulbright Foundation.

Simulation results presented in this paper has benefited from the Experimental
Infrastructure for Exploration of Exascale Computing (eX3), which is
financially supported by the Research Council of Norway under contract 270053.

## Author contributions statement

J.G.L. wrote the original code and conducted preliminary simulations. E.A.F.
conducted the simulations shown in this work, wrote the tests, analyzed data,
and prepared figures. J.S.D. and H.F. contributed to software development and
conducted numerical and performance testing. C.T.L. assisted in the original
development of SMART and worked on mesh processing. M.E.R. and P.R. designed
and supervised the project and wrote sections of the manuscript. All authors
reviewed the manuscript.

## Competing interests

The authors declare no competing interests.

##  Appendix A A general framework for multi-domain and multi-species reaction
and transport

SMART is designed to represent and solve coupled systems of nonlinear ordinary
and partial differential equations, describing reactions and/or transport due
to diffusion, convection or drift in multi-domain geometries including
subdomains of co-dimension zero or one (volume-surface problems). This section
defines the range and scope of SMART by describing its foundational
mathematical modelling framework. SMART uses a model specification framework
inspired by the Systems Biology Markup Language[29], consisting of
compartments, species, reactions, and parameters. Throughout this mathematical
overview, we refer to objects within SMART that correspond to different terms
in the equations.

###  A.1 Notation

  * •

$\mathcal{M}$: index/label set for domains $\Omega^{m}$ ($m\in\mathcal{M}$).

  * •

$\mathcal{Q}$: index/label set for the surfaces $\Gamma^{q}$
($q\in\mathcal{Q}$).

  * •

$\mathcal{I}$, $\mathcal{I}^{m}$, $\mathcal{I}^{q}$: index/label sets for the
species, the species in $\Omega^{m}$ and species on $\Gamma^{q}$,
respectively.

  * •

$\mathcal{K}$, $\mathcal{K}^{m}$,
$\mathcal{K}^{q},\mathcal{K}^{mq},\mathcal{K}^{mqn}$: index/label sets for the
reactions, reactions in $\Omega^{m}$, reactions on $\Gamma^{q}$, volume-
surface reactions between species in $\Omega^{m}$ and $\Gamma^{q}$, and
volume-surface-volume reactions between species in $\Omega^{m}$, $\Gamma^{q}$,
and $\Omega^{n}$.

###  A.2 Multi-domain and multi-surface geometry representation

We consider an open topologically $D$-dimensional manifold
$\Omega\subset\mathbb{R}^{d}$ for $D\leqslant d=1,2,3$, and assume that
$\Omega$ is partitioned into $|\mathcal{M}|$ open and disjoint _domains_
$\Omega^{m}\subset\mathbb{R}^{d}$:

$\Omega=\bigcup_{m\in\mathcal{M}}\Omega^{m},$ (S1)

each with (internal or external) boundary $\partial\Omega^{m}$, and boundary
$\partial\Omega$. In the SMART framework, each domain $\Omega^{m}$ is referred
to as a volume compartment. The _boundary normal_ $\mathbf{n}^{m}$ on
$\partial\Omega^{m}$ is the outward pointing normal vector field to the
boundary. The interface $\Gamma^{mn}$ between domains $\Omega^{m}$ and
$\Omega^{n}$ is defined as the intersection of the closure of the domains:

$\Gamma^{mn}=\overline{\Omega^{m}}\cap\overline{\Omega^{n}},$ (S2)

for $m,n\in\mathcal{M}$, and may or may not be empty. Each interface
$\Gamma^{mn}$ may be further partitioned into one or more _surfaces_ :

$\Gamma^{mn}=\bigcup_{q\in\mathcal{Q}^{mn}}\Gamma^{q},$ (S3)

where each $\Gamma^{q}$ is referred to as a surface compartment within SMART.
We denote the total set of surfaces by
$\mathcal{Q}=\cup_{m,n\in\mathcal{M}}\cup\mathcal{Q}_{0}$ where
$\mathcal{Q}_{0}$ also includes (sub)surfaces on the boundary of $\Omega$. For
each surface, we define the surface-to-neighbors map from the surface index
$q$ to its neighboring domain indices $m,n$: $q\mapsto m,n$. If the surface is
a real boundary surface, the map returns the single neighbor $q\mapsto m,$.

Moreover, let time $t\in[0,T]$ for $T>0$.

###  A.3 Coupled multi-domain reaction-transport equations

Different physical processes may occur on domains or on surfaces. For
instance, in the context of computational neuroscience at the cellular level,
different species and processes dominate in the dendrites and axons, and on
the plasma membrane versus at post-synaptic densities. In general, SMART is
designed to model different species coexisting, interacting and moving within
a (sub)domain or (sub)surface, between (sub)domains and across (sub)surfaces.
Within the SMART framework, information about each species is stored within an
object, and the equations describing their interactions are stored within
reaction objects. Any other quantities involved in the reaction equations
(e.g., reaction rates and binding affinities) are defined separately and
stored within parameter objects.

#### Domain species

We define and denote the set of species coexisting in $\Omega^{m}$ by
$\mathcal{I}^{m}$ for each (sub)domain $\Omega^{m}$. For each $m$, the species
$i\in\mathcal{I}^{m}$ with concentrations $u_{i}^{m}=u_{i}^{m}(x,t)$ for
$x\in\Omega^{m}$ and $t\in(0,T]$ satisfy reaction-transport equations of the
form

$\partial_{t}u_{i}^{m}+\mathcal{T}_{i}^{m}(u_{i}^{m})-f_{i}^{m}(u^{m})=0\qquad\text{
in }\Omega^{m},$ (S4)

where $\partial_{t}$ is the time-derivative, $\mathcal{T}_{i}^{m}$ defines
transport terms, and $f_{i}^{m}$ are volume reactions i.e. reactions within
the domain, typically given by non-linear and non-trivial relations involving
one or more of the species $u^{m}=\\{u_{i}^{m}\\}_{i\in\mathcal{I}^{m}}$. In
the case of transport by diffusion alone,

$\mathcal{T}_{i}^{m}(u)=-\nabla\cdot(D_{i}^{m}\nabla(u)),$ (S5)

where $\nabla\cdot$ is the spatial divergence operator, $\nabla$ is the
spatial gradient, and $D_{i}^{m}$ is the diffusion coefficient of species $i$
in domain $\Omega^{m}$ which may be heterogeneous i.e., spatially varying
and/or anisotropic i.e., tensor-valued. Remaining relative to $\Omega^{m}$, on
(all) surfaces $\Gamma^{q}\subseteq\Gamma^{mn}$, we assume that the flux of
species $i$ is governed by a relation of the form:

$D_{i}\nabla u_{i}^{m}\cdot n^{m}-R_{i}^{q}(u^{m},u^{n},v^{q})=0\qquad\text{
on }\Gamma^{q},$ (S6)

where $R_{i}^{q}$ defines surface fluxes, and the surface concentrations
$v^{q}$ are defined below.

#### Surface species

The surface concentrations $v^{q}=\\{v_{j}^{q}\\}_{j\in\mathcal{I}_{q}}$,
entering in (S6) above, are defined on $\Gamma^{q}\subseteq\Gamma^{mn}$,
either as prescribed fields or via surface equations as follows: find
$v^{q}=v^{q}(x,t)$ for $x\in\Gamma^{q}$, $t\in(0,T]$ such that

$\partial_{t}v_{j}^{q}+\mathcal{T}_{i}^{q}(v_{i}^{q})-g_{j}^{q}(u^{m},u^{n},v^{q})=0\qquad\text{
on }\Gamma^{q},$ (S7)

where $g_{j}^{q}$ are surface reactions for each species
$j\in\mathcal{I}^{q}$. In the case of transport by surface diffusion, with
surface diffusion coefficient $D_{j}^{q}$ for species $j$, surface gradient
$\nabla_{S}$ and surface divergence $\nabla_{S}\cdot$, we have:

$\mathcal{T}_{j}^{q}(v)=-\nabla_{S}\cdot(D_{j}^{q}\nabla_{S}v).$ (S8)

#### Boundary conditions

For any (sub)boundary $\Gamma^{q}\subset\partial\Omega^{m}$ such that
$\Gamma^{q}\subseteq\partial\Omega$, the following boundary conditions are
prescribed:

$D_{i}^{m}\nabla u_{i}^{m}\cdot n^{m}-R_{i}^{q}(u^{m},v^{q})=0,$ (S9)

which can be seen as a special case of Equation S6 in which there is no
adjacent volume compartment. For any species $v^{q}$ on any surface
$\Gamma^{q}$ with non-empty boundary, zero flux boundary conditions are
prescribed:

$D_{j}^{q}\nabla_{S}v_{j}^{q}\cdot n^{q}=0,$ (S10)

where $n^{q}$ is the normal to the boundary line.

#### Initial conditions

Initial conditions are required for any $u_{i}^{m}$ for $m\in\mathcal{M}$,
$i\in\mathcal{I}^{m}$, and $v_{i}^{q}$ for $q\in\mathcal{Q}$,
$i\in\mathcal{I}^{q}$.

###  A.4 SMART model generation from reaction specifications

Here, we outline the conventions used internally by SMART to convert reaction
specifications into associated terms in PDEs and boundary conditions. In
SMART, each reaction has several “flux" objects associated with it - one for
every reactant and every product. The manner in which these fluxes are
computed depends on the type of reaction, generally classified as “volume",
“surface", or “volume-surface" (Figure 2).

Considering the set of all expressions associated with volume reactions within
compartment $\Omega^{m}$,
$\mathcal{F}^{m}=\\{\mathcal{F}^{m}_{k}\\}_{k\in\mathcal{K}^{m}}$, the flux
associated with a reactant or product species $i$ for reaction $k$ is given
by:

$f_{i,k}^{m}=\mathcal{S}_{i,k}^{m}\mathcal{F}^{m}_{k},$ (S11)

where $\mathcal{S}_{i,k}^{m}$ is the stoichiometric coefficient, which is
positive for products and negative for reactants and whose magnitude matches
the number of molecules generated or consumed.

Surface reactions are handled analogously. The set of all surface reactions
within compartment $\Gamma^{q}$ is
$\mathcal{G}^{q}=\\{\mathcal{G}^{q}_{k}\\}_{k\in\mathcal{K}^{q}}$, and the
flux associated with a reactant or product species $j$ for reaction $k$ is
given by:

$g_{j,k}^{q}=\mathcal{S}_{j,k}^{q}\mathcal{G}^{q}_{k},$ (S12)

where $\mathcal{S}_{i,k}^{m}$ is the stoichiometric coefficient as before.

Volume-surface reactions generally contribute to both PDEs and boundary
conditions as follows. Considering the set of all expressions associated with
volume-surface reaction $k$ between compartments $\Omega^{m}$ and
$\Gamma^{q}$,
$\mathcal{R}^{mq}=\\{\mathcal{R}^{mq}_{k}\\}_{k\in\mathcal{K}^{mq}}$, the flux
associated with a reactant or product volume species $i$ for reaction $k$ is
given by:

$R_{i,k}^{mq}=\alpha_{i,k}^{mq}\mathcal{S}_{i,k}^{mq}\mathcal{R}^{mq}_{k},$
(S13)

where $\mathcal{S}_{i,k}^{m}$ is the stoichiometric coefficient and
$\alpha_{i,k}^{mq}$ is a scaling factor converting molecular flux (molecules
per unit surface per unit time) into the dimensionally equivalent units of
volume concentration multiplied by length per unit time. For any surface
species $j$ that acts as a product or reactant, such reactions have an
associated flux:

$g_{j,k}^{mq}=\mathcal{S}_{j,k}^{mq}\mathcal{R}^{mq}_{k},$ (S14)

where $\mathcal{S}_{j,k}^{mq}$ is the stoichiometric coefficient.

The remaining case of a surface-mediated reaction between two adjacent volume
compartments (“volume-surface-volume" reaction) is a natural extension of the
above “volume-surface" case. The set of expressions for all such reactions
between species in $\Omega^{m}$, $\Gamma^{q}$, and $\Omega^{n}$ is
$\mathcal{R}^{mqn}=\\{\mathcal{R}^{mqn}_{k}\\}_{k\in\mathcal{K}^{mqn}}$ and
the contribution to reactants and products are given by replacing all cases of
$mq$ with $mqn$ in Equations S13 and S14.

All of these fluxes are included in Equations S4, S6, S7 and S9 by summing the
associated contributions for a given species. In particular, the total
reaction term $f^{m}_{i}$ in Equation S4 is given by summing the contribution
from all volume reactions involving species $i$, $\mathcal{K}_{i}^{m}$:

$f^{m}_{i}=\sum_{k\in\mathcal{K}_{i}^{m}}{f_{i,k}^{m}}.$ (S15)

The total reaction term $R^{q}_{i}$ in Equation S6 or Equation S9 is given by
summing the contribution from all volume-surface or volume-surface-volume
reactions involving species $i$ ($\mathcal{K}_{i}^{mq}$,
$\mathcal{K}_{i}^{mqn}$):

$R^{q}_{i}=\sum_{k\in\mathcal{K}_{i}^{mq}}{R_{i,k}^{mq}}+\sum_{k\in\mathcal{K}_{i}^{mqn}}{R_{i,k}^{mqn}},$
(S16)

where the second term is equal to zero when applied to Equation S9, where
there is no adjacent volume compartment.

Finally, the total reaction term $g^{q}_{j}$ in Equation S7 is given by
summing the contribution from all surface reactions involving species $j$,
$\mathcal{K}_{j}^{q}$ as well as any volume-surface or volume-surface-volume
reactions involving species $j$ ($\mathcal{K}_{j}^{mq}$,
$\mathcal{K}_{j}^{mqn}$):

$g^{q}_{j}=\sum_{k\in\mathcal{K}_{j}^{q}}{g_{j,k}^{q}}+\sum_{k\in\mathcal{K}_{j}^{mq}}{g_{j,k}^{mq}}+\sum_{k\in\mathcal{K}_{j}^{mqn}}{g_{j,k}^{mqn}}.$
(S17)

####  A.4.1 Example of assembly from reactions

As a simple case of model generation in SMART, we consider the bulk-surface
reaction from Rangamani et al. [42], which is also included as Example 2 in
the SMART documentation. In brief, this model involves a volume species $A$ in
the cytosolic domain $\Omega_{cyto}$ that can bind to the surface species $X$
on the plasma membrane domain $\Gamma_{PM}$ to form surface species $B$
(Figure S1A). Treating $A$ and $X$ as reactants and $B$ as a product, the rate
of this volume-surface reaction is defined according to mass conservation:

$\mathcal{R}_{1}^{cyto-PM}=k_{on}u_{A}v_{X}-k_{off}v_{B}$ (S18)

From Equation S13, the boundary flux for species $A$ is then:

$R_{A,1}^{cyto-PM}=-\alpha_{A,1}^{cyto-PM}\mathcal{R}_{1}^{cyto-PM},$ (S19)

where $\alpha_{A,1}^{cyto-PM}$ depends on the chosen units. For instance, in a
3D/2D system with volume concentration units of $\upmu$M and surface
concentration units of molecules/$\upmu$m 2, $\alpha_{A,1}^{cyto-
PM}=\frac{1\times 10^{21}}{N_{A}}$ , where $N_{A}$ is Avogadro’s number. From
Equation S14, the surface reaction rate for $X$ is:

$g_{X,1}^{cyto-PM}=-\mathcal{R}_{1}^{cyto-PM},$ (S20)

and the surface reaction rate for $B$ is:

$g_{B,1}^{cyto-PM}=\mathcal{R}_{1}^{cyto-PM}.$ (S21)

In terms of these expressions, the system of PDEs is:

$\displaystyle\frac{\partial{u_{A}}}{\partial{t}}$
$\displaystyle=D_{A}\nabla^{2}u_{A}\quad\text{in}\;\Omega_{cyto}$ (S22)
$\displaystyle D_{A}(\textbf{n}\cdot\nabla u_{A})$
$\displaystyle=R_{A,1}^{cyto-PM}\quad\text{on}\;\Gamma_{PM}$ (S23)
$\displaystyle\frac{\partial{v_{X}}}{\partial{t}}$
$\displaystyle=D_{X}\nabla^{2}v_{X}+g_{X,1}^{cyto-
PM}\quad\text{on}\;\Gamma_{PM}$ (S24)
$\displaystyle\frac{\partial{v_{B}}}{\partial{t}}$
$\displaystyle=D_{B}\nabla^{2}v_{B}+g_{B,1}^{cyto-
PM}\quad\text{on}\;\Gamma_{PM}.$ (S25)

Figure S1: Simple example of a bulk-surface reaction in a cell. A) Schematic
summarizing the reaction of volume species A with surface species X to form
surface species B. B) Normalized concentration over time for each of the three
species.

##  Appendix B Numerical approximation and computational strategy

SMART solves the multi-domain reaction-transport equations outlined in Section
A via finite difference discretizations in time and a finite element
discretization in space, using a _monolithic_ approach to address the system
couplings. We describe this monolithic approach here, considering linear
diffusion-only transport relations for concreteness and to illustrate linear-
nonlinear strategies.

###  B.1 Monolithic solution of the multi-domain reaction-diffusion equations

Formally, we consider the Sobolev spaces $H^{1}(\Omega^{m})$,
$m\in\mathcal{M}$ of square-integrable functions on $\Omega^{m}$ with square-
integrable weak derivatives, and analogously for $H^{1}(\Gamma^{q})$, $q\in
Q$, as well as the vector function spaces
$H^{1}(\Omega^{m},\mathbb{R}^{d})=H^{1}(\Omega)^{d}$. For solving the multi-
domain reaction-diffusion equations, we introduce two product spaces $U$ and
$V$ consisting of bulk fields and surface fields, respectively:

$U=\bigotimes_{m\in\mathcal{M}}H^{1}(\Omega^{m};\mathbb{R}^{|\mathcal{I}^{m}|}),\quad
V=\bigotimes_{q\in\mathcal{Q}}H^{1}(\Gamma^{q};\mathbb{R}^{|\mathcal{I}^{q}|}).$
(S26)

To represent the solution fields $u,v$ comprised of the separate bulk and
surface components for the different species, we label

$u=\\{u^{m}\\}_{m\in\mathcal{M}}=\\{\\{u_{i}^{m}\\}_{i\in\mathcal{I}^{m}}\\}_{m\in\mathcal{M}},\quad
v=\\{v^{q}\\}_{q\in\mathcal{Q}}=\\{\\{v_{i}^{q}\\}_{i\in\mathcal{I}^{q}}\\}_{q\in\mathcal{Q}}.$
(S27)

By standard techniques (integrating by test functions and integrating by
parts), we rephrase the multi-domain reaction-diffusion equations in
variational form, resulting in the following coupled variational problem: find
$u\in U$ and $v\in V$ such that for all $\phi\in U$ and $\psi\in V$:

$F(u,v;\phi)+G(u,v;\psi)=0,$ (S28)

where both forms $F$ and $G$ are composed of sums over domains or surfaces and
species:

$F(u,v;\phi)=\sum_{m\in\mathcal{M}}\sum_{i\in\mathcal{I}^{m}}F_{i}^{m}(u,v;\phi_{i}^{m}),\qquad
G(u,v;\psi)=\sum_{q\in\mathcal{Q}}\sum_{i\in\mathcal{I}^{q}}G_{i}^{q}(u,v;\psi_{i}^{q}).$
(S29)

Furthermore, with the $L^{2}(O)$-inner product over any given domain
$O\subset\Omega$ defined as

$\langle a,b\rangle_{O}=\int_{O}a\cdot b\,\mathrm{d}x,$

we have defined

$F_{i}^{m}(u,v,\phi_{i}^{m})=\langle\partial_{t}u_{i}^{m},\phi_{i}^{m}\rangle_{\Omega^{m}}+\langle
D_{i}^{m}\nabla u_{i}^{m},\nabla\phi_{i}^{m}\rangle_{\Omega^{m}}-\langle
f_{i}^{m}(u^{m}),\phi_{i}^{m}\rangle_{\Omega^{m}}-\sum_{q\in\mathcal{Q}^{mn}}\langle
R_{i}^{q}(u^{m},u^{n},v^{q}),\phi_{i}^{m}\rangle_{\Gamma^{q}}$ (S30)

and, for any $q\in\mathcal{Q}^{mn}$ for given $\Omega^{m}$ and $\Omega^{n}$
interfacing $\Omega^{m}$ via $\Gamma^{q}$, finally:

$G_{i}^{q}(u,v,\psi_{i}^{q})=\langle\partial_{t}v_{i}^{q},\psi_{i}^{q}\rangle_{\Gamma^{q}}+\langle
D_{i}^{q}\nabla_{S}v_{i}^{q},\nabla_{S}\psi_{i}^{q}\rangle_{\Gamma^{q}}-\langle
g_{i}^{q}(u^{m},u^{n},v^{q}),\psi_{i}^{q}\rangle_{\Gamma^{q}}.$ (S31)

###  B.2 Discretization in time and space

SMART discretizes in time via an implicit (first-order) Euler scheme with time
steps $0=t_{0}<t_{1}<\ldots<t_{n-1}<t_{n}=T$ with timestep
$\tau_{n}=t_{n}-t_{n-1}$. To discretize in space, we employ a conforming
finite element mesh $\mathcal{T}$ defined relative to the parcellation of
$\Omega$ into domains and surfaces, such that
$\mathcal{T}^{m}\subseteq\mathcal{T}$ defines a submesh of $\Omega^{m}$ for
$m\in\mathcal{M}$, and $\mathcal{G}^{q}$ defines a (co-dimension one) mesh of
$\Gamma^{q}$ for $q\in\mathcal{Q}$. We define the finite element spaces of
continuous piecewise linears $\mathcal{P}_{1}$ defined relative to each
(sub)mesh, and define the product spaces

$U_{h}=\bigotimes_{m\in\mathcal{M}}\mathcal{P}_{1}^{|\mathcal{I}_{m}|}(\mathcal{T}^{m}),\qquad
V_{h}=\bigotimes_{q\in\mathcal{Q}}\mathcal{P}_{1}^{|\mathcal{I}_{q}|}(\mathcal{G}^{q}).$
(S32)

For each time step $t_{i}$, given discrete solutions $u_{h}^{-}$ and
$v_{h}^{-}$ at the previous time $t_{i-1}$, we then solve the nonlinear,
coupled time-discrete problem: find $u_{h}\in U_{h}$ and $v_{h}\in V_{h}$

$H(u_{h},v_{h})=F_{\tau_{n}}(u_{h},v_{h},\phi)+G_{\tau_{n}}(u_{h},v_{h},\psi)=0$
(S33)

for all $\phi\in U_{h}$, $\psi\in V_{h}$, where $F_{\tau_{n}}$ and
$G_{\tau_{n}}$ are defined by the forms $F$ and $G$ after implicit Euler time-
discretization and depend on the given $u_{h}^{-}$ and $v_{h}^{-}$.

###  B.3 Nonlinear solution algorithms

By default, the monolithic nonlinear discrete system (S33) is solved by
Newton-Raphson iteration with a symbolically derived discrete Jacobian.

##  Appendix C Software implementation

The SMART abstractions and algorithms are implemented via the open source and
generally available FEniCS finite element software package (2022-version)
[23]. FEniCS supports high-level specification of variational forms via the
Unified Form Language (UFL) [59], symbolic differentiation of variational
forms e.g. for derivation for Jacobians, automated assembly of general and
nonlinear forms over finite element meshes, and high-performance linear
algebra and nonlinear solvers via e.g. PETSc [31]. Below, we describe the
abstractions used by SMART to translate user-level specifications into the
variational form of the reaction-transport equations summarized above.

###  C.1 SMART model components

A model in SMART consists of a collection of four classes of objects -
compartments, species, reactions, and parameters.

A compartment is initialized by a name, dimensionality, associated length
units, and marker value in the parent mesh. Each compartment is linked to a
parent mesh by marker values stored in discrete functions defined over mesh
entities, $M_{facet}$ defined over all mesh facets (faces for 2D parent mesh,
edges for 3D parent mesh) and $M_{cell}$ defined over all mesh cells. For
instance, in a 3D/2D model, $M_{cell}$ stores a marker value associated with
each tetrahedron in the geometry; if the cytosol is linked to marker value
“1", then the cytosol is comprised of all tetrahedra with value “1" in
$M_{cell}$. All information about the compartment topology is implicit in the
mesh functions, but evaluation can be sped up by users explicitly specifying
which compartments are nonadjacent.

A species is initialized by a name, initial condition, units, diffusion
coefficient, diffusion coefficient units, and compartment. Importantly, each
species is associated with only one compartment; a single molecule that moves
between compartments will be specified by two species (e.g., cytosolic calcium
and ER calcium). The initial condition can be specified as a constant value or
as a string giving a spatially dependent expression. The unknown symbols in an
expression can include spatial coordinates $x,y,z$ or, in the case of surface
species, the curvature $curv$.

A reaction is initialized by a name, a list of reactants, a list of products,
equations for the forward and reverse reaction rates, and Python dictionaries
mapping strings from the equations to parameters and species initialized
separately. The reactant and product lists also specify the stoichiometry of
each reaction, e.g. , if a product is listed twice, it has a stoichiometric
coefficient of 2. Reactions are integrated into the variational forms of each
PDE as detailed above, and SMART carries out internal checks to ensure
consistency in units using the Python package pint [60]. If no forward or
reverse reaction strings are specified, SMART defaults to assuming mass action
kinetics, with an on-rate parameter on and an off-rate parameter off. For non-
mass-action kinetics, SMART supports a range of non-algebraic functions,
including trigonometric functions, exponential and logarithmic functions, and
absolute value and sign functions. Furthermore, reactions can be defined as
curvature sensitive by including $curv$ in the reaction expression.

Finally, a parameter is initialized by a name, value, and units. The value of
a parameter can be either a scalar or a string specifying a time and/or space-
dependent expression, where the free symbols should only include $x,y,z,t$. A
time-dependent parameter can also be loaded from a text file as shown in
Example 5 of the SMART documentation. If a time-dependent parameter with no
dependence on species concentrations is used to define a flux, the user may
optionally set `pre_integration` to `True`. In this case, SMART uses a pre-
integrated expression or a numerical approximation to update parameter values
while ensuring consistent mass transfer; that is, given a time-dependent flux
$f(t)$, flux values are updated as follows:

$f(t+\tau_{n})=\tau_{n}^{-1}\int_{t}^{t+\tau_{n}}\,f(T)dT$ (S34)

When using this updating scheme together with implicit Euler time integration,
the amount of mass transferred is independent of the time step.

###  C.2 Solver specifications and adaptive time-stepping

SMART uses PETSc4py for all linear and nonlinear solves [61]. Linear systems
are solved iteratively using Krylov solvers with field-split biconjugate
gradient preconditioning. The nonlinear system is solved using the Newton line
search SNES solver within PETSc. Default relative tolerances in each case are
set to $1\times 10^{-5}$, but these and other solver parameters can be readily
adjusted as illustrated in Example 6 of the SMART documentation.

By default, SMART takes uniform time-steps matching the initial specification.
However, there is an optional function for adaptive time-stepping using the
number of Newton iterations required to converge at the previous step. The
rules were set ad hoc and are summarized in Table S1. The reported “time-step
adjustment factor" $\zeta$ is the multiplicative factor used to alter the
time-step at the end of a given time step, i.e., the new time-step
$\tau_{new}$ is given by a in terms of the previous time-step $\tau_{old}$ as
$\tau_{new}=\zeta\tau_{old}$.

Table S1: Adaptive time stepping rules in SMART. Number of | Time-step
---|---
Nonlinear Iterations | Adjustment Factor
0-1 | 1.1
2-4 | 1.05
5-10 | 1.0
11-20 | 0.8
>20 | 0.5

###  C.3 Axisymmetric models

SMART offers the ability to assume axisymmetry given a 2D mesh. Given a 2D
mesh in the $r-z$ plane, the symmetry axis is assumed to be $r=0$. When this
feature is turned on in SMART, the variational forms are adjusted accordingly,
with the only necessary change being the multiplication by $r$ in each
integral (factors of $2\pi$ cancel throughout):

$\displaystyle F_{i}^{m}(u,v,\phi_{i}^{m})$ $\displaystyle=\langle
r\partial_{t}u_{i}^{m},\phi_{i}^{m}\rangle_{\Omega^{m}}+\langle
rD_{i}^{m}\nabla u_{i}^{m},\nabla\phi_{i}^{m}\rangle_{\Omega^{m}}-\langle
rf_{i}^{m}(u^{m}),\phi_{i}^{m}\rangle_{\Omega^{m}}-\sum_{q\in\mathcal{Q}^{mn}}\langle
rR_{i}^{q}(u^{m},u^{n},v^{q}),\phi_{i}^{m}\rangle_{\Gamma^{q}}$ (S35)
$\displaystyle G_{i}^{q}(u,v,\psi_{i}^{q})$ $\displaystyle=\langle
r\partial_{t}v^{q},\psi_{i}^{q}\rangle_{\Gamma^{q}}+\langle
rD_{i}^{q}\nabla_{S}v_{i}^{q},\nabla_{S}\psi_{i}^{q}\rangle_{\Gamma^{q}}-\langle
rg_{i}^{q}(u^{m},u^{n},v^{q}),\psi_{i}^{q}\rangle_{\Gamma^{q}}.$ (S36)

In this case, the compartments must be initialized as a 2D/1D system, as shown
in the mechanotransduction code and Example 3 with axisymmetry in the SMART
repository.

##  Appendix D Numerical testing

###  D.1 Numerical testing of protein phosphorylation model

As an idealized test case for SMART, we considered the simple system described
by Meyers et al., in which a single protein is phosphorylated at the plasma
membrane with rate $k_{kin}$ and dephosphorylated through the cytosol with
rate $k_{p}$ [41]. The SMART specifications are summarized in LABEL:tab:phos-
compartments, LABEL:tab:phos-species, LABEL:tab:phos-reactions and
LABEL:tab:phos-params.

Table S2: Compartments in phosphorylation model. Mesh statistics are given for the coarsest mesh ($h=3.0$) | Dimensionality | Species | Vertices | Cells | Marker value | Size
---|---|---|---|---|---|---
PM | 2 | 0 | 242 | 400 | 10 | $8.000\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
Cyto | 3 | 1 | 363 | 1200 | 1 | $8.000\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
Table S3: Species in phosphorylation model. | Compartment | $D(\unit{\micro\squared\per})$ | Initial condition
---|---|---|---
Aphos | Cyto | 10 | $1.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
Table S4: Reactions in phosphorylation model. | Reactants | Products | Equation | Type
---|---|---|---|---
r1 | [] | [‘Aphos’] | $(VolSA)k_{kin}(A_{tot}-A_{phos})$ | volume_surface
r2 | [‘Aphos’] | [] | $k_{dephos}A_{phos}$ | volume
Table S5: Parameters in phosphorylation model. | Value/Equation | Description
---|---|---
A_tot | $1.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic A
k_kin | $5.000\text{\times}{10}^{01}\text{\,}{\mathrm{s}}^{-1}$ | Rate constant for A phosphorylation
VolSA | $5.000\text{\times}{10}^{-02}\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | Cytosolic volume to PM surface area ratio
k_p | $1.000\text{\times}{10}^{01}\text{\,}{\mathrm{s}}^{-1}$ | Rate constant for A dephosphorylation

Given the two reactions summarized above, we can define the following flux
terms in accordance with Section A.4:

$\displaystyle R_{A_{phos},r1}^{cyto-PM}$
$\displaystyle=\alpha_{A_{phos},r1}^{cyto-
PM}[VolSA]k_{kin}(A_{tot}-u_{A_{phos}}),$ (S37) $\displaystyle
f_{A_{phos},r2}^{cyto}$ $\displaystyle=-k_{p}A_{phos}.$ (S38)

Accordingly the governing equation and boundary conditions are:

$\displaystyle\frac{\partial{u_{A_{phos}}}}{\partial{t}}$
$\displaystyle=D_{A_{phos}}\nabla^{2}u_{A_{phos}}+f_{A_{phos},r2}^{cyto}\quad\text{in}\;\Omega_{Cyto},$
(S39) $\displaystyle\quad D_{A_{phos}}(\textbf{n}\cdot\nabla u_{A_{phos}})$
$\displaystyle=R_{A_{phos},r1}^{cyto-PM}\quad\text{on}\;\Gamma_{PM}.$ (S40)

Due to the simplicity of this problem, a closed form solution exists for the
steady state concentration in 1D. This 1D solution approximates the 3D case in
which a cell is modeled as a very thin sheet, where the solution is used along
the thin cell dimension in the $z$ direction. A cell of thickness $\Delta z$
with one region of membrane at $z=0$ and the other at $z=\Delta z$ is then
predicted to have the following concentration profile at steady state:

$\displaystyle u_{A_{phos},SS}$
$\displaystyle=A_{tot}C_{1}\left(\exp\left(\frac{\Delta
z-z}{z_{0}}\right)+\exp\left(\frac{z}{z_{0}}\right)\right),$ (S41)
$\displaystyle\text{where}\qquad C_{1}$
$\displaystyle=\frac{k_{phos}z_{0}}{D_{A_{phos}}\left(\exp\left(\frac{\Delta
z}{z_{0}}\right)-1\right)+k_{phos}z_{0}\left(1+\exp\left(\frac{\Delta
z}{z_{0}}\right)\right)}$ (S42) $\displaystyle\text{and}\qquad z_{0}$
$\displaystyle=\sqrt{\frac{D_{A_{phos}}}{k_{dephos}}}.$ (S43)

We ran these simulations for a thin slab where the thickness is ten times
smaller than the length in the $x$ and $y$ directions, testing the effects of
time-step refinement, mesh refinement, and changes in the diffusion
coefficient. As an error metric, we computed the $\mathcal{L}_{2}$ norm of the
difference between the analytical steady state solution and our numerical
solution at $t=1$ s ($u_{A_{phos},num})$:

$\left\lvert\left\lvert
u_{A_{phos},SS}-u_{A_{phos},num}\right\rvert\right\rvert_{\mathcal{L}2}=\sqrt{\leavevmode\nobreak\
\iiint\limits_{\Omega_{Cyto}}{(u_{A_{phos},SS}-u_{A_{phos},num})^{2}d\Omega}}$
(S44)

Table S6: $\mathcal{L}_{2}$ error (in $\upmu$M) between computed and analytical solution for different time steps and mesh sizes for $D=10.0\leavevmode\nobreak\ \unit{\micro\squared\per}$. |  $h=0.375$ $\upmu$m |  $h=0.75$ $\upmu$m |  $h=1.5$ $\upmu$m |  $h=3$ $\upmu$m
---|---|---|---|---
$\tau=0.64$ s | $8.87\text{\times}{10}^{-02}\text{\,}$ | $8.58\text{\times}{10}^{-02}\text{\,}$ | $1.06\text{\times}{10}^{-01}\text{\,}$ | $3.68\text{\times}{10}^{-01}\text{\,}$
$\tau=0.32$ s | $1.97\text{\times}{10}^{-02}\text{\,}$ | $2.56\text{\times}{10}^{-02}\text{\,}$ | $8.77\text{\times}{10}^{-02}\text{\,}$ | $3.96\text{\times}{10}^{-01}\text{\,}$
$\tau=0.16$ s | $5.58\text{\times}{10}^{-03}\text{\,}$ | $2.24\text{\times}{10}^{-02}\text{\,}$ | $9.14\text{\times}{10}^{-02}\text{\,}$ | $4.04\text{\times}{10}^{-01}\text{\,}$
$\tau=0.08$ s | $5.69\text{\times}{10}^{-03}\text{\,}$ | $2.31\text{\times}{10}^{-02}\text{\,}$ | $9.21\text{\times}{10}^{-02}\text{\,}$ | $4.06\text{\times}{10}^{-01}\text{\,}$
$\tau=0.04$ s | $5.78\text{\times}{10}^{-03}\text{\,}$ | $2.32\text{\times}{10}^{-02}\text{\,}$ | $9.22\text{\times}{10}^{-02}\text{\,}$ | $4.06\text{\times}{10}^{-01}\text{\,}$
$\tau=0.02$ s | $5.80\text{\times}{10}^{-03}\text{\,}$ | $2.32\text{\times}{10}^{-02}\text{\,}$ | $9.23\text{\times}{10}^{-02}\text{\,}$ | $4.06\text{\times}{10}^{-01}\text{\,}$
$\tau=0.01$ s | $5.80\text{\times}{10}^{-03}\text{\,}$ | $2.32\text{\times}{10}^{-02}\text{\,}$ | $9.23\text{\times}{10}^{-02}\text{\,}$ | $4.06\text{\times}{10}^{-01}\text{\,}$
Table S7: $\mathcal{L}_{2}$ error (in $\upmu$M) between computed and analytical solution for different time steps and mesh sizes for $D=100.0\leavevmode\nobreak\ \unit{\micro\squared\per}$. |  $h=0.375$ $\upmu$m |  $h=0.75$ $\upmu$m |  $h=1.5$ $\upmu$m |  $h=3$ $\upmu$m
---|---|---|---|---
$\tau=0.64$ s | $9.64\text{\times}{10}^{-02}\text{\,}$ | $9.58\text{\times}{10}^{-02}\text{\,}$ | $9.35\text{\times}{10}^{-02}\text{\,}$ | $8.41\text{\times}{10}^{-02}\text{\,}$
$\tau=0.32$ s | $2.09\text{\times}{10}^{-02}\text{\,}$ | $2.04\text{\times}{10}^{-02}\text{\,}$ | $1.98\text{\times}{10}^{-02}\text{\,}$ | $4.27\text{\times}{10}^{-02}\text{\,}$
$\tau=0.16$ s | $2.53\text{\times}{10}^{-03}\text{\,}$ | $2.93\text{\times}{10}^{-03}\text{\,}$ | $9.33\text{\times}{10}^{-03}\text{\,}$ | $4.78\text{\times}{10}^{-02}\text{\,}$
$\tau=0.08$ s | $5.89\text{\times}{10}^{-04}\text{\,}$ | $2.40\text{\times}{10}^{-03}\text{\,}$ | $9.86\text{\times}{10}^{-03}\text{\,}$ | $4.89\text{\times}{10}^{-02}\text{\,}$
$\tau=0.04$ s | $6.06\text{\times}{10}^{-04}\text{\,}$ | $2.48\text{\times}{10}^{-03}\text{\,}$ | $9.96\text{\times}{10}^{-03}\text{\,}$ | $4.90\text{\times}{10}^{-02}\text{\,}$
$\tau=0.02$ s | $6.19\text{\times}{10}^{-04}\text{\,}$ | $2.49\text{\times}{10}^{-03}\text{\,}$ | $9.97\text{\times}{10}^{-03}\text{\,}$ | $4.90\text{\times}{10}^{-02}\text{\,}$
$\tau=0.01$ s | $6.22\text{\times}{10}^{-04}\text{\,}$ | $2.49\text{\times}{10}^{-03}\text{\,}$ | $9.98\text{\times}{10}^{-03}\text{\,}$ | $4.90\text{\times}{10}^{-02}\text{\,}$
Table S8: $\mathcal{L}_{2}$ error (in $\upmu$M) between computed and analytical solution for different time steps and mesh sizes for $D=1000.0\leavevmode\nobreak\ \unit{\micro\squared\per}$. |  $h=0.375$ $\upmu$m |  $h=0.75$ $\upmu$m |  $h=1.5$ $\upmu$m |  $h=3$ $\upmu$m
---|---|---|---|---
$\tau=0.64$ s | $9.72\text{\times}{10}^{-02}\text{\,}$ | $9.71\text{\times}{10}^{-02}\text{\,}$ | $9.69\text{\times}{10}^{-02}\text{\,}$ | $9.49\text{\times}{10}^{-02}\text{\,}$
$\tau=0.32$ s | $2.11\text{\times}{10}^{-02}\text{\,}$ | $2.10\text{\times}{10}^{-02}\text{\,}$ | $2.08\text{\times}{10}^{-02}\text{\,}$ | $1.92\text{\times}{10}^{-02}\text{\,}$
$\tau=0.16$ s | $2.65\text{\times}{10}^{-03}\text{\,}$ | $2.59\text{\times}{10}^{-03}\text{\,}$ | $2.48\text{\times}{10}^{-03}\text{\,}$ | $4.40\text{\times}{10}^{-03}\text{\,}$
$\tau=0.08$ s | $3.04\text{\times}{10}^{-04}\text{\,}$ | $3.27\text{\times}{10}^{-04}\text{\,}$ | $9.39\text{\times}{10}^{-04}\text{\,}$ | $4.90\text{\times}{10}^{-03}\text{\,}$
$\tau=0.04$ s | $6.54\text{\times}{10}^{-05}\text{\,}$ | $2.38\text{\times}{10}^{-04}\text{\,}$ | $9.91\text{\times}{10}^{-04}\text{\,}$ | $5.02\text{\times}{10}^{-03}\text{\,}$
$\tau=0.02$ s | $5.95\text{\times}{10}^{-05}\text{\,}$ | $2.48\text{\times}{10}^{-04}\text{\,}$ | $1.00\text{\times}{10}^{-03}\text{\,}$ | $5.04\text{\times}{10}^{-03}\text{\,}$
$\tau=0.01$ s | $6.12\text{\times}{10}^{-05}\text{\,}$ | $2.50\text{\times}{10}^{-04}\text{\,}$ | $1.01\text{\times}{10}^{-03}\text{\,}$ | $5.04\text{\times}{10}^{-03}\text{\,}$

###  D.2 Extra details on the numerical testing of biological test cases

In our numerical testing of the mechanotransduction and dendritic spine
examples, we considered refined versions of the reference mesh. For the model
of mechanotransduction, we tested refined versions of the 2D mesh used to
model an axisymmetric spread cell on a circular contact region. The mesh
statistics and total degrees of freedom (DOFs) are given in Table S9.

Table S9: Refined meshes for numerical testing of the mechanotransduction model. Refinements | Vertices | Triangles | DOFs
---|---|---|---
0 | 308 | 548 | 3078
1 | 1163 | 2192 | 11058
2 | 4517 | 8768 | 41772
3 | 17801 | 35072 | 162216

In the case of the dendritic spine example, we tested a range of different
mesh refinements and time steps, as summarized in Figure S2. Overall, we found
the effects of time-step refinement were almost identical across three
different mesh refinements. We observe almost identical solutions for time
steps 0.001 s or less, justifying our choice of 0.001 s time steps for the
results shown in Figure 4. At the smallest time step tested, the error due to
mesh refinement is very small, as indicated by the minimal differences across
mesh refinements. As shown in Figure 6E in the main text, mesh refinement does
influence the spatial solution, but this effect appears to be quite small at
the starting resolution of the dendritic spine mesh.

Figure S2: Summary of calcium dynamics in dendritic spine upon mesh refinement
and time-step refinement. Curves are plotted for three different mesh
refinements and 5 different time steps, as indicated. The calcium peak for the
coarsest mesh is indicated by the horizontal dashed line throughout, showing
the reduction in calcium peak upon mesh refinement.

The mesh statistics for each refinement are provided in Table S10.

Table S10: Refined meshes for numerical testing of the dendritic spine model. Refinements | Vertices | Tetrahedra | DOFs
---|---|---|---
0 | 18649 | 106979 | 49194
1 | 146024 | 855832 | 323328
2 | 1154871 | 6846656 | 2299090

##  Appendix E Mesh generation for mechanotransduction model

Cell geometries on micropatterned substrates were generated in Gmsh [57]. As a
reference geometry, we used the axisymmetric cell shape previously defined in
VCell [33], where the plasma membrane surface satisfies

$\frac{\left(1-\frac{z^{4}}{2000+z^{4}}\right)(r^{2}+z^{2})+0.4(r^{2}+(z+9.72)^{2})z^{4}}{15+z^{4}}=169$
(S45)

and the nuclear membrane surface satisfies

$\left(\frac{r}{5.3}\right)^{2}+\left(\frac{z-4.8}{2.4}\right)^{2}=1,$ (S46)

where $r$ and $z$ are cylindrical coordinates. These implicit boundary
representations were constructed by solving for $r$ at each $z$ value using
solveset in Sympy [62].

Shapes with rectangular or star-shaped contact regions were constructed by
introducing an angular dependence to the reference shape. In particular, we
scale the $r$ value at the surface as a function of the cylindrical coordinate
$\theta$. Considering an implicitly (as above) or explicitly defined contour
as a function of arc length $r(s),z(s)$, the scaled contour at a given value
of $\theta$ is $r(s)T(\theta),z(s)$. We additionally require that this shape
conserves cell volume, resulting in the constraint:

$\int_{0}^{2\pi}{T^{2}(\theta)d\theta}=2\pi.$ (S47)

For the circular contact region, $T(\theta)=1$ and the above constraint is
automatically satisfied.

For a rectangular contact region with sides of length $a$ and $b$:

$T_{rect}(\theta)=\begin{cases}\frac{a}{2\cos{\theta}}&-\arctan(\frac{b}{a})<\theta\leq\arctan(\frac{b}{a})\\\
\frac{b}{2\sin{\theta}}&\arctan(\frac{b}{a})<\theta\leq\pi-\arctan(\frac{b}{a})\\\
-\frac{a}{2\cos{\theta}}&\pi-\arctan(\frac{b}{a})<\theta\leq\pi+\arctan(\frac{b}{a})\\\
-\frac{b}{2\sin{\theta}}&\pi+\arctan(\frac{b}{a})<\theta\leq
2\pi-\arctan(\frac{b}{a})\\\ \end{cases},$ (S48)

resulting in the constraint $ab=\pi$. For our example, we chose
$a=\sqrt{0.6\pi}$ and $b=\sqrt{\frac{\pi}{0.6}}$ for an aspect ratio of 0.6.

For a star-shaped contact region:

$T_{star}(\theta)=T_{0}+T_{1}\cos(5\theta),$ (S49)

resulting in the constraint $2T_{0}^{2}+T_{1}^{2}=2$. For our example, we
chose $T_{0}=0.98$ and $T_{1}=0.2814$.

Using these shapes directly results in sharp corners and unrealistic
artifacts, so we applied additional smoothing steps without altering the cell
volume. We first constructed an additional shape with an elliptical contact
region; $T(\theta)$ for an ellipse with axes $c,d$ is:

$T_{ellipse}(\theta)=\frac{cd}{\sqrt{\left(d\cos(\theta)\right)^{2}+\left(c\sin(\theta)\right)^{2}}},$
(S50)

and the associated constraint is $cd=1$. We choose the aspect ratio of the
ellipse to be $\frac{\text{range}(x_{contact})}{\text{range}(y_{contact})}$,
where $x_{contact}$ and $y_{contact}$ are the set of coordinates within the
star or rectangular contact region. Accordingly, this shape is a smoother
option with a contact region that very roughly approximates the desired
contour. To smooth out the cell shape away from the substrate, we can define
new $T$ functions that also depend on $z$:

$\displaystyle T_{rect,smooth}^{2}(\theta,z)$
$\displaystyle=\frac{z-z_{max}}{z_{max}}T_{rect}^{2}+\frac{z}{z_{max}}T_{ellipse}^{2},$
(S51) $\displaystyle T_{star,smooth}^{2}(\theta,z)$
$\displaystyle=\frac{z-z_{max}}{z_{max}}T_{star}^{2}+\frac{z}{z_{max}}T_{ellipse}^{2},$
(S52)

where $z_{max}$ is the height of the cell contour. Using these definitions,
the desired contour shape at the surface is perfectly preserved and cell
volume is conserved. In the case of the rectangular mesh, we also rounded off
the sharp corners of the contact region, replacing them with arcs of radius
0.2 $\upmu$m. The nucleus shape remained unaltered across all simulations.

##  Appendix F Model specifications for all examples and biological test cases

Here, we summarize each model in terms of its compartments, species, reactions
and parameters. We note that these tables were directly output from SMART
using the print_to_latex functions associated with each class of container,
with some minor modifications for readability.

###  F.1 Mechanotransduction model

Table S11: Compartments in mechanotransduction model. | Dimensionality | Species | Vertices | Cells | Marker value | Size
---|---|---|---|---|---|---
Circular contact region |  |  |  |  |  | ($\int{2\pi rdr}$)
Cyto | 2 | 11 | 7920 | 15360 | 1 | $1.925\text{\times}{10}^{03}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
PM | 1 | 2 | 305 | 304 | 10 | $1.295\text{\times}{10}^{03}\text{\,}\mathrm{\SIUnitSymbolMicro m}$
Nuc | 2 | 1 | 649 | 1159 | 2 | $2.823\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
NM | 1 | 2 | 125 | 124 | 12 | $2.345\text{\times}{10}^{02}\text{\,}\mathrm{\SIUnitSymbolMicro m}$
Rectangular contact region |  |  |  |  |  | (1/4 symmetry)
Cyto | 3 | 11 | 17479 | 84954 | 1 | $4.794\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
PM | 2 | 2 | 5113 | 10030 | 10 | $3.314\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
Nuc | 3 | 1 | 1668 | 6406 | 2 | $7.044\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
NM | 2 | 2 | 901 | 1716 | 12 | $5.856\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
Star-shaped contact region |  |  |  |  |  | (1/10 symmetry)
Cyto | 3 | 11 | 7687 | 34466 | 1 | $1.899\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
PM | 2 | 2 | 2158 | 4110 | 10 | $1.360\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
Nuc | 3 | 1 | 799 | 2759 | 2 | $2.818\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
NM | 2 | 2 | 380 | 674 | 12 | $2.342\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
Table S12: Species in mechanotransduction model. | Compartment |  $D$ ($\unit{\micro\squared\per}$) | Initial Condition
---|---|---|---
Emod∗ | PM | 0 |  $70\left(\dfrac{1-\text{sign}(z-0.0001)}{2}\right)$ GPa
pFAK | Cyto | 10 | $3.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
RhoA_GDP | Cyto | 1 | $1.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
RhoA_GTP | PM | 0.3 | $3.360\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
ROCK_A | Cyto | 75 | $0.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
mDia_A | Cyto | 1 | $0.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
Myo_A | Cyto | 0.8 | $1.500\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
LIMK_A | Cyto | 10 | $1.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
Cofilin_NP | Cyto | 10 | $1.800\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
FActin | Cyto | 0.6 | $1.790\text{\times}{10}^{01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
GActin | Cyto | 13.37 | $4.824\text{\times}{10}^{02}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
LaminA | NM | 0.001 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
NPC_A | NM | 0.001 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
YAPTAZ | Cyto | 19 | $7.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
YAPTAZ_phos | Cyto | 19 | $2.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
YAPTAZ_nuc | Nuc | 19 | $7.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
∗Emod represents the substrate stiffness and also specifies the location of the substrate ($z=0$). |  |  |
Table S13: Reactions in mechanotransduction model. | Reactants | Products | Equation | Type
---|---|---|---|---
a1 | [] | [‘pFAK’] | $(cyto_{Convert})([FAK_{tot}]-[pFAK])(\dfrac{E_{mod}k_{sf}}{C+Emod}+k_{f})$ | vol-surf
a2 | [‘pFAK’] | [] | $k_{df}[pFAK]$ | vol
a3 | [‘RhoA_GDP’] | [‘RhoA_GTP’] | $[RhoA_{GDP}]cyto_{Convert}k_{fkrho}(\gamma_{const}pFAK^{5}+1)$ | vol-surf
|  |  | $-[RhoA_{GTP}]k_{drho}$ |
a4 | [] | [‘ROCK_A’] | $[RhoA_{GTP}]k_{rrho}([ROCK_{tot}]-[ROCK_{A}])$ | vol-surf
a5 | [‘ROCK_A’] | [] | $[ROCK_{A}]k_{drock}$ | vol
b1 | [] | [‘mDia_A’] | $[RhoA_{GTP}]k_{mrho}([mDia_{tot}]-[mDia_{A}])$ | vol-surf
b2 | [‘mDia_A’] | [] | $k_{dmdia}[mDia_{A}]$ | vol
b3 | [] | [‘Myo_A’] | $k_{mr}([Myo_{tot}]-[Myo_{A}])\cdot$ | vol
|  |  | $\left([ROCK_{A}]\epsilon\dfrac{1+\tanh\left(sc1([ROCK_{A}]-[ROCK_{B}])\right)}{2}+1\right)$ |
|  |  | $-[Myo_{A}]k_{dmy}$ |
b4 | [] | [‘LIMK_A’] | $k_{lr}([LIMK_{tot}]-[LIMK_{A}])\cdot$ | vol
|  |  | $\left([ROCK_{A}]\tau\dfrac{1+\tanh\left(sc1([ROCK_{A}]-[ROCK_{B}])\right)}{2}+1\right)$ |
|  |  | $-[LIMK_{A}]k_{dl}$ |
b5 | [] | [‘Cofilin_NP’] | $k_{turnover}([Cofilin_{tot}]-[Cofilin_{NP}])$ | vol
|  |  | $-\dfrac{[Cofilin_{NP}][LIMK_{A}]k_{catCof}}{Cofilin_{NP}+k_{mCof}}$ |
b6 | [‘GActin’] | [‘FActin’] | $[GActin]k_{ra}\cdot$ | vol
|  |  | $\left(\alpha[mDia_{A}]\dfrac{1+\tanh\left(sc1([mDia_{A}]-[mDia_{B}])\right)}{2}+1\right)$ |
|  |  | $-[FActin]([Cofilin_{NP}]k_{fc1}+k_{dep})$ |
c1 | [‘YAPTAZ_phos’] | [‘YAPTAZ’] | $[YAPTAZ_{phos}]([FActin][Myo_{A}]k_{CY}+k_{CN})$ | vol
|  |  | $-[YAPTAZ]k_{NC}$ |
c3 | [] | [‘LaminA’] | $\dfrac{[FActin]^{2.6}k_{fl}p([LaminA_{tot}]-[LaminA])}{C_{LaminA}+[FActin]^{2.6}p}-[LaminA]k_{rl}$ | vol-surf
c4 | [] | [‘NPC_A’] | $[FActin][LaminA][Myo_{A}]k_{fNPC}([NPC_{tot}]-[NPC_{A}])$ | vol-surf
|  |  | $-[NPC_{A}]k_{rNPC}$ |
c5 | [‘YAPTAZ’] | [‘YAPTAZ_nuc’] | $[YAPTAZ]([NPC_{A}]k_{in2}+k_{insolo})-[YAPTAZ_{nuc}]k_{out}$ | vol-surf-vol
Table S14: Parameters in mechanotransduction model. Values were adopted directly from those reported in [33] | Value | Description
---|---|---
FAK_tot | $1.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic FAK
k_f | $1.500\text{\times}{10}^{-02}\text{\,}{\mathrm{s}}^{-1}$ | rate constant for baseline FAK phosphorylation
k_sf | $3.790\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | substrate-stiffness-dependent FAK phosphorylation rate constant
C | $3.250\text{\times}{10}^{00}\text{\,}\mathrm{kPa}$ | critical substrate stiffness
cytoConvert | $1.825\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | vol/SA ratio
k_df | $3.500\text{\times}{10}^{-02}\text{\,}{\mathrm{s}}^{-1}$ | FAK dephosphorylation rate constant
k_fkrho | $1.680\text{\times}{10}^{-02}\text{\,}{\mathrm{s}}^{-1}$ | RhoA activation rate constant
gammaConst | $7.756\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-5}$ | Scaling factor
k_drho | $6.250\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | RhoA deactivation rate constant
k_rrho | $6.480\text{\times}{10}^{-01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | ROCK activation rate constant
ROCK_tot | $1.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic ROCK
k_drock | $8.000\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | ROCK deactivation rate constant
k_mrho | $2.000\text{\times}{10}^{-03}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | mDia activation rate constant
mDia_tot | $8.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic mDia
k_dmdia | $5.000\text{\times}{10}^{-03}\text{\,}{\mathrm{s}}^{-1}$ | mDia deactivation rate constant
Myo_tot | $5.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic myosin
k_mr | $3.000\text{\times}{10}^{-02}\text{\,}{\mathrm{s}}^{-1}$ | Myosin activation rate constant
ROCK_B | $3.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Critical activated ROCK concentration
epsilon | $3.600\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}$ | Scaling factor
sc1 | $2.000\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}$ | Scaling factor
k_dmy | $6.700\text{\times}{10}^{-02}\text{\,}{\mathrm{s}}^{-1}$ | Myosin deactivation rate constant
LIMK_tot | $2.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic LIMK
k_lr | $7.000\text{\times}{10}^{-02}\text{\,}{\mathrm{s}}^{-1}$ | LIMK activation rate constant
tau | $5.549\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}$ | Scaling factor
k_dl | $2.000\text{\times}{10}^{00}\text{\,}{\mathrm{s}}^{-1}$ | LIMK deactivation rate constant
Cofilin_tot | $2.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic cofilin
k_turnover | $4.000\text{\times}{10}^{-02}\text{\,}{\mathrm{s}}^{-1}$ | Cofilin dephosphorylation rate constant
k_catCof | $3.400\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | Cofilin phosphorylation rate constant
k_mCof | $4.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Critical cofilin concentration
k_ra | $4.000\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | Actin polymerization rate constant
alpha | $5.000\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}$ | Scaling factor
mDia_B | $1.650\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Critical activated mDia concentration
k_dep | $3.500\text{\times}{10}^{00}\text{\,}{\mathrm{s}}^{-1}$ | Actin depolymerization rate constant
k_fc1 | $4.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | Cofilin-mediated actin polymerization rate constant
k_CN | $5.600\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | Baseline YAP/TAZ dephosphorylation rate constant
k_CY | $7.600\text{\times}{10}^{-04}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ | stress-fiber-mediated YAP/TAZ dephosphorylation rate constant
k_NC | $1.400\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | YAP/TAZ phosphorylation rate constant
LaminA_tot | $3.500\text{\times}{10}^{03}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$ | Total lamin A in NM
k_fl | $4.600\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | Lamin A dephosphorylation rate constant
p | $9.000\text{\times}{10}^{-06}\text{\,}\mathrm{kPa}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-2.600}$ | Scaling factor for cell stiffness
C_LaminA | $1.000\text{\times}{10}^{02}\text{\,}\mathrm{kPa}$ | Critical cell stiffness for lamin A dephosphorylation
k_rl | $1.000\text{\times}{10}^{-03}\text{\,}{\mathrm{s}}^{-1}$ | Lamin A phosphorylation rate constant
NPC_tot | $6.500\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$ | Total NPC density in NM
k_fNPC | $2.800\text{\times}{10}^{-07}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ | NPC opening rate constant
k_rNPC | $8.700\text{\times}{10}^{00}\text{\,}{\mathrm{s}}^{-1}$ | NPC closing rate constant
k_insolo | $1.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | NPC-independent YAP/TAZ nuclear import rate
k_in2 | $1.000\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | NPC-dependent YAP/TAZ nuclear import rate
k_out | $1.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | YAP/TAZ nuclear export rate

###  F.2 Dendritic spine calcium model

Table S15: Compartments in dendritic spine model. | Dimensionality | Species | Vertices | Cells | Marker value | Size
---|---|---|---|---|---|---
Cyto | 3 | 2 | 17742 | 81206 | 1 | $6.463\text{\times}{10}^{-01}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
PM | 2 | 3 | 1661 | 3243 | 10 | $6.236\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
SA | 3 | 1 | 8727 | 25773 | 2 | $2.692\text{\times}{10}^{-02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
SAm | 2 | 0 | 7820 | 15672 | 12 | $2.387\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
Table S16: Species in dendritic spine model. | Compartment |  $D$ ($\unit{\micro\squared\per}$) | Initial Condition
---|---|---|---
Ca | Cyto | 220 | $1.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
NMDAR | PM | 0 | $$1.000\text{\times}{10}^{00}\text{\,}$^{*}$
VSCC | PM | 0 | $\dfrac{1+\text{sign}(z+0.25)}{2}$
Bf | PM | 0 | $7.957\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
Bm | Cyto | 20 | $2.000\text{\times}{10}^{01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
CaSA | SA | 45 | $6.000\text{\times}{10}^{01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
∗NMDAR restricted to postsynaptic density. |  |  |
Table S17: Reactions in dendritic spine model. | Reactants | Products | Equation | Type
---|---|---|---|---
a1 | [] | [‘Ca^2+’] | $G_{NMDAR}(J0_{NMDAR})[NMDAR](V_{m}-V_{rev})$ | volume_surface
a2 | [] | [‘Ca^2+’] | $J_{VSCC}[VSCC]$ | volume_surface
a3 | [‘Ca^2+’] | [] | $\dfrac{100n_{PMr}\left(\dfrac{[Ca^{2+}]^{5}(Vmax_{hr23})}{[Ca^{2+}]^{5}+(Km_{hr23})^{5}}+\dfrac{[Ca^{2+}]^{2}(Vmax_{lr23})}{[Ca^{2+}]^{2}+(Km_{lr23})^{2}}\right)}{1+\dfrac{(Kme)(Prtote)}{([Ca^{2+}]+(Kme))^{2}}+\dfrac{(Kmx)(Prtotex)}{([Ca^{2+}]+(Kmx))^{2}}}$ | volume_surface
a4 | [‘Ca^2+’] | [] | $\dfrac{1000n_{PMr}\dfrac{(Vmax_{r22})[Ca^{2+}]}{[Ca^{2+}]+(Km_{r22})}}{1+\dfrac{(Kme)(Prtote)}{([Ca^{2+}]+(Kme))^{2}}+\dfrac{(Kmx)(Prtotex)}{([Ca^{2+}]+(Kmx))^{2}}}$ | volume_surface
a5 | [‘Ca^2+’, ‘Bf’] | [] | $(kBf_{on})[Bf][Ca^{2+}]-(kBf_{off})((Bf_{tot})-[Bf])$ | volume_surface
b1 | [‘Ca^2+’, ‘Bm’] | [] | $(kBm_{on})[Bm][Ca^{2+}]-(kBm_{off})((Bm_{tot})-[Bm])$ | volume
c1 | [‘Ca^2+’] | [‘CaSA’] | $\dfrac{1000n_{SAr}\dfrac{[Ca^{2+}]^{2}(Vmax_{r19})}{[Ca^{2+}]^{2}+(KP_{r19})^{2}}}{1+\dfrac{(Kme)(Prtote)}{([Ca^{2+}]+(Kme))^{2}}+\dfrac{(Kmx)(Prtotex)}{([Ca^{2+}]+(Kmx))^{2}}}$ | volume_surface_volume
c2 | [‘CaSA’] | [‘Ca^2+’] | $k_{leak}n_{SAr}([CaSA]-[Ca^{2+}])$ | volume_surface_volume
Table S18: Parameters in dendritic spine model. | Value/Equation | Description
---|---|---
n_PMr | $1.011\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | Experimental vol to surface area ratio
Vm | $\big{[}\left(55.51\exp\left(-\dfrac{t}{0.003}\right)+10.29\exp(-40.0t)\right)\dfrac{\text{sign}(t)+1}{2}$ | PM voltage
|  $+25\left(\exp\left(-\dfrac{t}{0.05}\right)-\exp\left(-\dfrac{t}{0.005}\right)\right)\dfrac{\text{sign}(t)+1}{2}-65\big{]}$ mV |
Vrev | $9.000\text{\times}{10}^{01}\text{\,}\mathrm{mV}$ | Reversal voltage for NMDAR
G_NMDAR | $1.377\text{\times}{10}^{05}\text{\,}{\mathrm{mV}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | NMDAR conductance
J0_NMDAR |  $-\dfrac{0.00456\left(\exp\left(-\dfrac{t}{0.05}\right)+\exp\left(-\dfrac{t}{0.05}\right)\right)(\text{sign}(t)+1)}{110.8\exp(-0.092V_{m})+1}$ $\upmu$m -2 |  Ca^2+ influx through NMDAR
J_VSCC | $(0.393-2.245\exp(-0.01044V_{m})(\exp(-34700t)-\exp(-3680t))\cdot$ |  Ca^2+ influx through VSCC
|  $(\text{sign}(t)+1)\dfrac{-1120(2V_{m}+1456)}{192970-85942\exp(0.01044Vm)}\text{s}^{-1}$ $\upmu$m -2 |
Prtote | $1.910\text{\times}{10}^{02}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ |  Ca^2+ buffer 1 concentration
Kme | $2.430\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ |  Ca^2+ buffer 1 affinity
Prtotex | $8.770\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ |  Ca^2+ buffer 2 concentration
Kmx | $1.390\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ |  Ca^2+ buffer 2 affinity
Vmax_lr23 | $1.130\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}\text{\,}{\mathrm{s}}^{-1}$ | PMCA conductance 1
Km_lr23 | $4.420\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | PMCA activation constant 1
Vmax_hr23 | $5.900\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}\text{\,}{\mathrm{s}}^{-1}$ | PMCA conductance 2
Km_hr23 | $4.420\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | PMCA activation constant 2
Vmax_r22 | $1.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}\text{\,}{\mathrm{s}}^{-1}$ | NCX conductance
Km_r22 | $1.000\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | NCX activation constant
kBf_on | $1.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | On-rate for fixed Ca^2+ buffer
kBf_off | $2.000\text{\times}{10}^{00}\text{\,}{\mathrm{s}}^{-1}$ | Off-rate for fixed Ca^2+ buffer
Bf_tot | $7.957\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total amount of fixed Ca^2+ buffer
kBm_on | $1.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | On-rate for mobile Ca^2+ buffer
kBm_off | $1.000\text{\times}{10}^{00}\text{\,}{\mathrm{s}}^{-1}$ | Off-rate for mobile Ca^2+ buffer
Bm_tot | $2.000\text{\times}{10}^{01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total conc. of mobile Ca^2+ buffer
n_SAr | $1.130\text{\times}{10}^{-02}\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | Cytosolic volume to SA surf. area ratio
Vmax_r19 | $1.140\text{\times}{10}^{02}\text{\,}\mathrm{\SIUnitSymbolMicro M}\text{\,}{\mathrm{s}}^{-1}$ | SERCA conductance
KP_r19 | $2.000\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | SERCA activation constant
k_leak | $1.608\text{\times}{10}^{-01}\text{\,}{\mathrm{s}}^{-1}$ | SR Ca^2+ leak rate constant

###  F.3 Calcium release unit model

Table S19: Compartments in CRU model. | Dimensionality | Species | Vertices | Cells | Marker value | Size
---|---|---|---|---|---|---
Cyto | 3 | 4 | 45557 | 204894 | 1 | $8.029\text{\times}{10}^{08}\text{\,}{\mathrm{nm}}^{3}$
SR | 3 | 2 | 12178 | 42030 | 2 | $3.438\text{\times}{10}^{06}\text{\,}{\mathrm{nm}}^{3}$
TTM | 2 | 0 | 13282 | 26552 | 10 | $7.190\text{\times}{10}^{05}\text{\,}{\mathrm{nm}}^{2}$
SRM | 2 | 1 | 9749 | 19518 | 12 | $4.129\text{\times}{10}^{05}\text{\,}{\mathrm{nm}}^{2}$
Table S20: Species in CRU model. | Compartment |  $D$ ($\unit{\micro\squared\per}$) | Initial Condition
---|---|---|---
Ca | Cyto | 220 | $1.400\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
ATP | Cyto | 140 | $4.547\text{\times}{10}^{02}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
CMDN | Cyto | 25 | $2.353\text{\times}{10}^{01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
TRPN | Cyto | 0 |  $\dfrac{56.8}{1+\exp\left(-\dfrac{(x+125)^{2}+(z+75)^{2}}{100^{2}}\right)}$ $\upmu$M ∗
CaSR | SR | 220 | $1.300\text{\times}{10}^{03}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
CSQN | SR | 25 | $2.110\text{\times}{10}^{03}\text{\,}\mathrm{\SIUnitSymbolMicro M}$
RyR | SRM | 0 | $$1.000\text{\times}{10}^{00}\text{\,}$^{**}$
∗ TRPN distribution chosen to approximate that in Ref. 12, using the distance from a central point in the T-tubule geometry rather than the distance to the junctional T-tubule boundary. |  |  |
∗∗ As in Ref. 12, RyR was localized to the SR-T-tubule junction; here defined as any region of the SR membrane less than 20 nm from the T-tubule membrane. |  |  |
Table S21: Reactions in CRU model. | Reactants | Products | Equation | Type
---|---|---|---|---
s1 | [‘CaSR’] | [‘Ca^2+’] | $(NO_{RyR})[RyR]g_{RyR}([CaSR]-[Ca^{2+}])$ | vol-surf-vol
s2∗ | [‘Ca^2+’] | [‘CaSR’] | $2(SR_{VtoA})S_{SERCA}\rho_{SERCA}\cdot$ | vol-surf-vol
|  |  | $\dfrac{[Ca^{2+}]^{2}\alpha_{1,plus}\alpha_{2,plus}-[CaSR]^{2}\alpha_{1,minus}\alpha_{2,minus}}{[Ca^{2+}]^{2}\alpha_{1,plus}+[CaSR]^{2}\alpha_{1,minus}+\alpha_{2,minus}+\alpha_{2,plus}}$ |
t1 | [] | [‘Ca^2+’] | $(Jmax)\left(1+(VFactor2)(ksat)\right)\cdot$ | vol-surf
|  |  | $\dfrac{(CaTT)(Nai)^{3}(VFactor1)-[Ca^{2+}](Nao)^{3}(VFactor2)}{[Ca^{2+}](denom2)+(denom1)}$ |
t2 | [‘Ca^2+’] | [] | $(Jp_{max})\dfrac{(Km_{pCa})^{2}}{[Ca^{2+}]^{2}+(Km_{pCa})^{2}}$ | vol-surf
t3 | [‘Ca^2+’] | [] | $(JCab)\left((Voltage)-0.5\dfrac{(R)(T)}{F}\log\left(\dfrac{(CaTT)}{[Ca^{2+}]}\right)\right)$ | vol-surf
b1 | [‘Ca^2+’, ‘ATP’] | [] | $(kon_{ATP})[ATP][Ca^{2+}]$ | vol
|  |  | $-(koff_{ATP})((ATP_{tot})-[ATP])$ |
b2 | [‘Ca^2+’, ‘CMDN’] | [] | $(kon_{CMDN})[CMDN][Ca^{2+}]$ | vol
|  |  | $-(koff_{CMDN})((CMDN_{tot})-[CMDN])$ |
b3 | [‘Ca^2+’, ‘TRPN’] | [] | $(kon_{TRPN})[Ca^{2+}][TRPN]$ | vol
|  |  | $-(koff_{TRPN})((TRPN_{tot})-[TRPN])$ |
b4 | [‘CaSR’, ‘CSQN’] | [] | $(kon_{CSQN})[CSQN][CaSR]$ | vol
|  |  | $-(koff_{CSQN})((CSQN_{tot})-[CSQN])$ |
∗In simulations without SERCA, this reaction was excluded from the model. |  |  |  |
Table S22: Parameters in CRU model. | Value | Description
---|---|---
ATP_tot | $4.550\text{\times}{10}^{02}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic ATP
CMDN_tot | $2.400\text{\times}{10}^{01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic CMDN
TRPN_tot | $7.000\text{\times}{10}^{01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic TRPN
CSQN_tot | $6.390\text{\times}{10}^{03}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Total cytosolic CSQN
CaTT | $1.800\text{\times}{10}^{03}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ |  Ca^2+ conc. in T-tubules
NO_RyR | $5.000\text{\times}{10}^{00}\text{\,}$ | Number of open RyR in release event
g_RyR | $2.620\text{\times}{10}^{04}\text{\,}\mathrm{nm}\text{\,}{\mathrm{s}}^{-1}$ | RyR conductance
S_SERCA | $1.500\text{\times}{10}^{00}\text{\,}$ | SERCA scale parameter
rho_SERCA | $7.500\text{\times}{10}^{01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Volume density of SERCA
SR_VtoA | $3.070\text{\times}{10}^{02}\text{\,}\mathrm{nm}$ | cytosolic volume to SA surf. area ratio 0
alpha1_plus | $1.067\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ | SERCA rate constant
alpha2_plus | $5.350\text{\times}{10}^{00}\text{\,}{\mathrm{s}}^{-1}$ | SERCA rate constant
alpha1_minus | $5.069\text{\times}{10}^{-05}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ | SERCA rate constant
alpha2_minus | $3.868\text{\times}{10}^{-02}\text{\,}{\mathrm{s}}^{-1}$ | SERCA rate constant
Voltage | $-8.200\text{\times}{10}^{-02}\text{\,}\mathrm{V}$ | TT voltage
R | $8.310\text{\times}{10}^{00}\text{\,}\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ | Universal gas constant
T | $2.9515\text{\times}{10}^{02}\text{\,}\mathrm{K}$ | Temperature
F | $9.650\text{\times}{10}^{04}\text{\,}\mathrm{C}\text{\,}{\mathrm{mol}}^{-1}$ | Faraday’s constant
Jmax | $4.083\text{\times}{10}^{02}\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}\mathrm{\SIUnitSymbolMicro M}\text{\,}{\mathrm{s}}^{-1}$ | Max. NCX flux
Nai | $1.420\text{\times}{10}^{04}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Cytosolic sodium conc.
Nao | $1.400\text{\times}{10}^{05}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ | Extracellular sodium conc.
VFactor1 | $3.233\text{\times}{10}^{-01}\text{\,}$ | NCX lumped parameter
VFactor2 | $8.142\text{\times}{10}^{00}\text{\,}$ | NCX lumped parameter
ksat | $2.700\text{\times}{10}^{-01}\text{\,}$ | NCX saturation parameter
denom1 | $3.852\text{\times}{10}^{16}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{4}$ | NCX lumped parameter
denom2 | $4.278\text{\times}{10}^{15}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{3}$ | NCX lumped parameter
Jp_max | $4.974\text{\times}{10}^{00}\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}\mathrm{\SIUnitSymbolMicro M}\text{\,}{\mathrm{s}}^{-1}$ | Max. Ca^2+ pump flux
Km_pCa | $2.890\text{\times}{10}^{-01}\text{\,}\mathrm{\SIUnitSymbolMicro M}$ |  Ca^2+ pump Ca^2+ affinity
JCab | $2.611\text{\times}{10}^{02}\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}\mathrm{\SIUnitSymbolMicro M}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{V}}^{-1}$ | TT Ca^2+ leak flux
kon_ATP | $2.250\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |  Ca^2+-ATP binding on rate
koff_ATP | $4.500\text{\times}{10}^{04}\text{\,}{\mathrm{s}}^{-1}$ |  Ca^2+-ATP binding off rate
kon_CMDN | $3.400\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |  Ca^2+-CMDN binding on rate
koff_CMDN | $2.380\text{\times}{10}^{02}\text{\,}{\mathrm{s}}^{-1}$ |  Ca^2+-CMDN binding off rate
kon_TRPN | $3.270\text{\times}{10}^{01}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |  Ca^2+-TRPN binding on rate
koff_TRPN | $1.960\text{\times}{10}^{01}\text{\,}{\mathrm{s}}^{-1}$ |  Ca^2+-TRPN binding off rate
kon_CSQN | $1.020\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro M}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |  Ca^2+-CSQN binding on rate
koff_CSQN | $6.500\text{\times}{10}^{04}\text{\,}{\mathrm{s}}^{-1}$ |  Ca^2+-CSQN binding off rate

###  F.4 ATP production in mitochondria

Table S23: Compartments in ATP generation model. | Dimensionality | Species | Vertices | Cells | Marker value | Size
---|---|---|---|---|---|---
IMS | 3 | 1 | 36667 | 128395 | 1 | $2.031\text{\times}{10}^{-02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
OM | 2 | 0 | 11514 | 23024 | 10 | $6.320\text{\times}{10}^{-01}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
Mat | 3 | 2 | 22974 | 85773 | 2 | $1.647\text{\times}{10}^{-02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{3}$
IM | 2 | 13 | 18708 | 37598 | 12 | $1.545\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$
∗In simulations with cristae-localized membrane species, the IM was split into two regions as shown in Figure 5B. The cristae portion had an area of 1.0197 $\upmu$m 2, with 28578 total surface elements. |  |  |  |  |  |
Table S24: Species in ATP generation model. | Compartment∗ |  $D$ ($\unit{\micro\squared\per}$) | Initial Condition ∗∗
---|---|---|---
E_Mat | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
E_IMS | IM | 0 | $$1.728\text{\times}{10}^{02}\text{\,}$/$2.618\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$$
E_Mat_H3Star | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
E_Mat_H3S | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
E_Mat_H3 | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
L | IM | 0 | $$1.066\text{\times}{10}^{04}\text{\,}$/$1.615\text{\times}{10}^{04}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$$
TL | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
LT | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
DL | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
LD | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
TLD | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
DLT | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
DLD † | IM | 0 | $0.000\text{\times}{10}^{00}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$
D_Mat | Mat | 15 | $7.200\text{\times}{10}^{-01}\text{\,}\mathrm{mM}$
T_Mat | Mat | 15 | $6.500\text{\times}{10}^{-01}\text{\,}\mathrm{mM}$
T_IMS | IMS | 15 | $3.250\text{\times}{10}^{-01}\text{\,}\mathrm{mM}$
∗In simulations with cristae-localized membrane species, all species listed as “IM" here were localized to the Cristae compartment instead. |  |  |
∗∗ Two values are given for some species - first value associated with uniform distribution of proteins in the IM, second value for cristae-localized proteins, ensuring consistency in total number of molecules. |  |  |
† Note that $TLT^{\prime}$ and $DLD^{\prime}$ introduced in Ref. 40 are not used here, as they are indistinguishable states that can be eliminated by substituting $TLT\leftarrow TLT+TLT^{\prime}$ and $DLD\leftarrow DLD+DLD^{\prime}$. |  |  |
Table S25: Reactions in ATP generation model. | Reactants | Products | Equation | Type
---|---|---|---|---
E1 | [‘E_Mat’] | [‘E_IMS’] | $k_{16}E_{Mat}-k_{61}E_{IMS}$ | surf
E2 | [‘E_IMS’] | [] | $E_{IMS}k_{65}-k_{56}E_{IMS,H3}\qquad^{*}$ | surf
E3 | [] | [‘E_Mat_H3Star’] | $k_{54}E_{IMS,H3}-E_{Mat,H3Star}k_{45}\qquad^{*}$ | surf
E4 | [] | [‘E_Mat_H3’] | $k_{52}E_{IMS,H3}-E_{Mat,H3}k_{25}\qquad^{*}$ | surf
E5 | [‘E_Mat_H3Star’, ‘D_Mat’] | [‘E_Mat_H3S’] | $k_{43}E_{Mat,H3Star}D_{Mat}-k_{34}E_{Mat,H3S}$ | vol-surf
E6 | [‘E_Mat_H3S’] | [‘E_Mat_H3’, ‘T_Mat’] | $k_{32}E_{Mat_{H}3S}-k_{23}E_{Mat,H3}T_{Mat}$ | vol-surf
E7 | [‘E_Mat_H3’] | [‘E_Mat’] | $k_{21}E_{Mat,H3}-k_{12}E_{Mat}$ | surf
L1 | [‘L’, ‘T_Mat’] | [‘LT’] | $k_{on,Tm}[L]T_{Mat}-k_{off,Tm}[LT]$ | vol-surf
L2 | [‘L’, ‘D_Mat’] | [‘LD’] | $k_{on,Dm}[L]D_{Mat}-k_{off,Dm}[LD]$ | vol-surf
L3 | [‘L’, ‘T_IMS’] | [‘TL’] | $k_{on,Ti}[L]T_{IMS}-k_{off,Ti}[TL]$ | vol-surf
L4 | [‘L’] | [‘DL’] | $D_{IMS}[L]k_{on,Di}-[DL]k_{off,Di}$ | surf
L5 | [‘TL’, ‘T_Mat’] | [] | $[TL]T_{Mat}k_{on,Tm}-k_{off,Tm}[TLT]\qquad^{**}$ | vol-surf
L6 | [‘TL’, ‘D_Mat’] | [‘TLD’] | $k_{on,Dm}[TL]D_{Mat}-k_{off,Dm}[TLD]$ | vol-surf
L7 | [‘LT’, ‘T_IMS’] | [] | $[LT]T_{IMS}k_{on,Ti}-k_{off,Ti}[TLT]\qquad^{**}$ | vol-surf
L8 | [‘LT’] | [‘DLT’] | $D_{IMS}[LT]k_{on,Di}-[DLT]k_{off,Di}$ | surf
L9 | [‘DL’, ‘T_Mat’] | [‘DLT’] | $k_{on,Tm}[DL]T_{Mat}-k_{off,Tm}[DLT]$ | vol-surf
L10 | [‘DL’, ‘D_Mat’] | [‘DLD’] | $k_{on,Dm}[DL]D_{Mat}-k_{off,Dm}[DLD]$ | vol-surf
L11 | [‘LD’, ‘T_IMS’] | [‘TLD’] | $k_{on,Ti}[LD]T_{IMS}-k_{off,Ti}[TLD]$ | vol-surf
L12 | [‘LD’] | [‘DLD’] | $D_{IMS}[LD]k_{on,Di}-[DLD]k_{off,Di}$ | surf
L13 | [‘DLT’] | [‘TLD’] | $k_{p}[DLT]-k_{cp}[TLD]$ | surf
V1 | [‘T_IMS’] | [‘T_cyto’] | $[VDAC]k_{vdac}(T_{IMS}-T_{cyto})$ | vol-surf
∗ Due to mass conservation, $E_{IMS,H3}$ is not treated as a separate variable, but solved for as $E_{IMS,H3}=E_{tot}-E_{Mat}-E_{IMS}-E_{Mat,H3Star}-E_{Mat,H3}-E_{Mat,H3S}$. |  |  |  |
∗∗ Due to mass conservation, $TLT$ is not treated as a separate variable, but solved for as $TLT=L_{tot}-L-LT-TL-LD-DL-DLD-DLT-TLD$ |  |  |  |
Table S26: Parameters for mitochondrial ATP generation. | Value | Description
---|---|---
E_tot | $1.73\text{\times}{10}^{02}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$ | Total ATP synthase
L_tot | $1.07\text{\times}{10}^{04}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$ | Total ANTs
D_IMS | $4.50\text{\times}{10}^{-02}\text{\,}\mathrm{mM}$ | IMS ADP concentration
k_16 | $1.48\text{\times}{10}^{05}\text{\,}{\mathrm{s}}^{-1}$ | ATP-synthase-associated rate constants
k_61 | $3.37\text{\times}{10}^{04}\text{\,}{\mathrm{s}}^{-1}$ |
k_65 | $3.97\text{\times}{10}^{03}\text{\,}{\mathrm{s}}^{-1}$ |
k_56 | $1.00\text{\times}{10}^{03}\text{\,}{\mathrm{s}}^{-1}$ |
k_54 | $1.00\text{\times}{10}^{02}\text{\,}{\mathrm{s}}^{-1}$ |
k_45 | $1.00\text{\times}{10}^{02}\text{\,}{\mathrm{s}}^{-1}$ |
k_52 | $1.00\text{\times}{10}^{-20}\text{\,}{\mathrm{s}}^{-1}$ |
k_25 | $5.85\text{\times}{10}^{-30}\text{\,}{\mathrm{s}}^{-1}$ |
k_43 | $2.00\text{\times}{10}^{03}\text{\,}{\mathrm{mM}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |
k_34 | $1.00\text{\times}{10}^{02}\text{\,}{\mathrm{s}}^{-1}$ |
k_32 | $5.00\text{\times}{10}^{03}\text{\,}{\mathrm{s}}^{-1}$ |
k_23 | $5.00\text{\times}{10}^{03}\text{\,}{\mathrm{mM}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |
k_21 | $4.00\text{\times}{10}^{01}\text{\,}{\mathrm{s}}^{-1}$ |
k_12 | $1.00\text{\times}{10}^{02}\text{\,}{\mathrm{s}}^{-1}$ |
koff_Tm | $4.00\text{\times}{10}^{04}\text{\,}{\mathrm{s}}^{-1}$ | ANT-associated rate constants
kon_Tm | $6.40\text{\times}{10}^{03}\text{\,}{\mathrm{mM}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |
koff_Ti | $2.00\text{\times}{10}^{02}\text{\,}{\mathrm{s}}^{-1}$ |
kon_Ti | $4.00\text{\times}{10}^{02}\text{\,}{\mathrm{mM}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |
koff_Dm | $4.00\text{\times}{10}^{04}\text{\,}{\mathrm{s}}^{-1}$ |
kon_Dm | $4.00\text{\times}{10}^{03}\text{\,}{\mathrm{mM}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |
koff_Di | $1.00\text{\times}{10}^{02}\text{\,}{\mathrm{s}}^{-1}$ |
kon_Di | $4.00\text{\times}{10}^{03}\text{\,}{\mathrm{mM}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ |
k_p | $9.20\text{\times}{10}^{01}\text{\,}{\mathrm{s}}^{-1}$ |
k_cp | $3.50\text{\times}{10}^{00}\text{\,}{\mathrm{s}}^{-1}$ |
k_vdac | $1.00\text{\times}{10}^{03}\text{\,}{\mathrm{mM}}^{-1}\text{\,}{\mathrm{s}}^{-1}$ | opening rate of VDACs
VDAC | $1.00\text{\times}{10}^{04}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-2}$ | Density of VDACs in OM

##  Appendix G Video captions

Video 1:  F-actin and YAP/TAZ dynamics in a cell on a circular micropattern.
F-actin concentration in the cytosol and YAP/TAZ concentration in the nucleus
plotted within a cell on a circular contact region with a radius of 13
$\upmu$m. Substrate stiffness matches that of a glass coverslip (70 GPa).
Video 2:  F-actin and YAP/TAZ dynamics in a cell on a rectangular
micropattern. F-actin concentration in the cytosol and YAP/TAZ concentration
in the nucleus plotted within a cell on a rectangular contact region with the
same contact area and cell volume as the cell in Video 1. Substrate stiffness
matches that of a glass coverslip (70 GPa).  Video 3:  F-actin and YAP/TAZ
dynamics in a cell on a star-shaped micropattern. F-actin concentration in the
cytosol and YAP/TAZ concentration in the nucleus plotted within a cell on a
star-shaped contact region with the same contact area and cell volume as the
cell in Video 1. Substrate stiffness matches that of a glass coverslip (70
GPa).  Video 4:  Calcium dynamics in a dendritic spine. Calcium influx into
the dendritic spine over the duration of a single voltage peak. Both cytosolic
and SR Ca^2+ concentrations are shown, with the SR in the spine head
accumulating Ca^2+ at later time points due to entry via SERCA.  Video 5:
Calcium dynamics in a cardiomyocyte CRU with SERCA. Cytosolic and SR calcium
concentrations during a single calcium release event from the SR. Ca^2+ is
slowly replenished in the SR after release due to entry through SERCA. A
smaller calcium flux occurs through the T-tubules (dark grey) and the two
mitochondria (light grey) act as a diffusional barrier.  Video 6:  Calcium
dynamics in a cardiomyocyte CRU without SERCA. Cytosolic and SR calcium
concentrations during a single calcium release event from the SR. In contrast
to Video 5, Ca^2+ remains depleted in the SR here as SERCA is not present in
the SR membrane.  Video 7:  ATP dynamics in a mitochondrion. ATP
concentrations in the matrix (innermost volume) and IMS (region between the
outer and inner membranes) are displayed over 100 ms of simulation. Rapid
initial decrease in both compartments is followed by a slower increase of IMS
ATP and a slow decrease in matrix ATP.

## References

  * [1] Yue, R. & Dutta, A.  Computational systems biology in disease modeling and control, review and perspectives.  _npj Systems Biology and Applications_ 8, 1–16, DOI: 10.1038/s41540-022-00247-4 (2022).
  * [2] Mogilner, A., Wollman, R. & Marshall, W. F.  Quantitative Modeling in Cell Biology: What Is It Good for?  _Developmental Cell_ 11, 279–287, DOI: 10.1016/j.devcel.2006.08.004 (2006).
  * [3] Kitano, H.  Computational systems biology.  _Nature_ 420, 206–210, DOI: 10.1038/nature01254 (2002).
  * [4] WAINWRIGHT, S. A.  Form and Function in Organisms.  _American Zoologist_ 28, 671–680, DOI: 10.1093/icb/28.2.671 (1988).
  * [5] Herman, M. A. _et al._ A Unifying Framework for Understanding Biological Structures and Functions Across Levels of Biological Organization.  _Integrative and Comparative Biology_ 61, 2038–2047, DOI: 10.1093/icb/icab167 (2022).
  * [6] Peddie, C. J. _et al._ Volume electron microscopy.  _Nature Reviews Methods Primers_ 2, 1–23, DOI: 10.1038/s43586-022-00131-9 (2022).