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

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 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 (   =  0 (   −   )) was assumed starting at  = 0.This realistic ER geometry was derived from electron micrographs of dendritic spines 30 .
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 (Ω  ) and 3 surface compartments (Γ  ), 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 Appendix 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    .Results are post-processed to examine changes in concentration over time and space.membrane surface area is increased for geometries with greater curvature at the contact regions (see Appendix E for details on generating these meshes).
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 to 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 µM, lower than the 8 µ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 µm 2 SA per µm 3 cytosol; original: 2.31 um 2 SA per um 3 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 µM near the T-tubule / SR junction and about 1 µ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.

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 5/39  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.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 .
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 Appendix 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 µm, 1.5 µm, 0.75 µm, and 0.375 µ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 , as indicated by the plateau in each curve (Figure 6B).This critical time-step is smaller for high values of , 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 L 2 norm of the error with respect to element size ℎ (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 µm 2 s −1 and treat all surface species as uniform and non-diffusing.We start with a time step of Δ = 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 Table S11).
Dendritic spine simulations were initially conducted on a mesh with ℎ  = 0.009 µm and ℎ  = 0.199 µm and a time step (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.
of Δ = 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 ℎ  = 0.00661 µm and ℎ  = 0.0383 µm, and our time step is adaptively adjusted from an initial value of Δ = 0.02 ms.When setting the diffusion coefficients for all volumetric nucleotide species to 150 µm 2 s −1 , 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).

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× increase in the number of cells for each refinement and a similar increase (7.8 − 7.9×) in the number of vertices.The corresponding numbers of degrees of freedom (, dimensions of the discrete solutions) are: 49 194, 323 328, 2 299 090; thus corresponding to a 6.6× and 7.1× increase between refinements.For each simulation, we use a 0.001 s timestep through  = 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 ), we observe that the simulation time scales close to log-linearly (between O () and O ( log )) (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).

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 9/39   S6, Table S7 and Table S8 for diffusion coefficients  10, 100, and 1000 µm 2 s −1 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.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 (Ω) and the lower dimensional boundaries or interfaces as the surfaces (Γ).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 Ω  ⊆ Ω are defined as a set of mesh cells with a common label or tag  ∈ M, and similarly exterior boundaries and interior interfaces Γ  are defined as a set of mesh facets again with a common tag  ∈ 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 Ω  and/or surface Γ  .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 Appendix 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 Ω  , and on or across the surfaces Γ  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 Appendix 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  𝑚  defined in each bulk subdomain Ω  and the concentration of each species    defined over each surface Γ  .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 13/39 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 -dependence in all integrals (Appendix 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  = 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 ( =   ), the ATP concentration   at 14/39 the next time point   +   was estimated as: where    is the ATP concentration in the IMS,   is the volume of cytosol immediately surrounding the mitochondrion,   is Avogadro's number, and other parameters are summarized in Appendix F. This updated value of   was then used in solving the PDE system, after which the estimated value was updated for consistency with the implicit Euler time discretization: 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 codimension 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.
• I, I  , I  : index/label sets for the species, the species in Ω  and species on Γ  , respectively.
• K, K  , K  , K  , K  : index/label sets for the reactions, reactions in Ω  , reactions on Γ  , volume-surface reactions between species in Ω  and Γ  , and volume-surface-volume reactions between species in Ω  , Γ  , and Ω  .

A.2 Multi-domain and multi-surface geometry representation
We consider an open topologically -dimensional manifold Ω ⊂ R  for  ⩽  = 1, 2, 3, and assume that each with (internal or external) boundary Ω  , and boundary Ω.In the SMART framework, each domain Ω  is referred to as a volume compartment.The boundary normal n  on Ω  is the outward pointing normal vector field to the boundary.The interface Γ  between domains Ω  and Ω  is defined as the intersection of the closure of the domains: for ,  ∈ M, and may or may not be empty.Each interface Γ  may be further partitioned into one or more surfaces: where each Γ  is referred to as a surface compartment within SMART.We denote the total set of surfaces by Q = ∪ ,∈ M ∪ Q 0 where Q 0 also includes (sub)surfaces on the boundary of Ω.For each surface, we define the surface-to-neighbors map from the surface index  to its neighboring domain indices , :  ↦ → , .If the surface is a real boundary surface, the map returns the single neighbor  ↦ → ,.Moreover, let time  ∈ [0,] for  > 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 Ω  by I  for each (sub)domain Ω  .For each , the species  ∈ I  with concentrations    =    (, ) for  ∈ Ω  and  ∈ (0,] satisfy reaction-transport equations of the form where   is the time-derivative, T   defines transport terms, and    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   = {   }  ∈ I  .In the case of transport by diffusion alone, where ∇• is the spatial divergence operator, ∇ is the spatial gradient, and    is the diffusion coefficient of species  in domain Ω  which may be heterogeneous i.e., spatially varying and/or anisotropic i.e., tensor-valued.Remaining relative to Ω  , on (all) surfaces Γ  ⊆ Γ  , we assume that the flux of species  is governed by a relation of the form: where    defines surface fluxes, and the surface concentrations   are defined below.

Surface species
The surface concentrations   = {   }  ∈ I  , entering in (S6) above, are defined on Γ  ⊆ Γ  , either as prescribed fields or via surface equations as follows: find   =   (, ) for  ∈ Γ  ,  ∈ (0,] such that where    are surface reactions for each species  ∈ I  .In the case of transport by surface diffusion, with surface diffusion coefficient    for species , surface gradient ∇  and surface divergence ∇  •, we have:

Boundary conditions
For any (sub)boundary Γ  ⊂ Ω  such that Γ  ⊆ Ω, the following boundary conditions are prescribed: which can be seen as a special case of Equation (S6) in which there is no adjacent volume compartment.For any species   on any surface Γ  with non-empty boundary, zero flux boundary conditions are prescribed: where   is the normal to the boundary line.

Initial conditions
Initial conditions are required for any    for  ∈ M,  ∈ I  , and    for  ∈ Q,  ∈ I  .

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 Ω  , F  = {F   }  ∈ K  , the flux associated with a reactant or product species  for reaction  is given by: where S  , 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 Γ  is G  = {G   }  ∈ K  , and the flux associated with a reactant or product species  for reaction  is given by: where S  , 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  between compartments Ω  and Γ  , R  = {R   }  ∈ K  , the flux associated with a reactant or product volume species  for reaction  is given by: where S  , is the stoichiometric coefficient and   , 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  that acts as a product or reactant, such reactions have an associated flux: where S  , 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 Ω  , Γ  , and Ω  is R  = {R   }  ∈ K  and the contribution to reactants and products are given by replacing all cases of  with  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  𝑚  in Equation (S4) is given by summing the contribution from all volume reactions involving species , K   : The total reaction term    in Equation (S6) or Equation (S9) is given by summing the contribution from all volume-surface or volume-surface-volume reactions involving species  (K   , K   ): 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    in Equation (S7) is given by summing the contribution from all surface reactions involving species , K   as well as any volume-surface or volume-surface-volume reactions involving species  (K   , K   ):

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  in the cytosolic domain Ω  that can bind to the surface species  on the plasma membrane domain Γ  to form surface species  (Figure S1A).
Treating  and  as reactants and  as a product, the rate of this volume-surface reaction is defined according to mass conservation: R From Equation (S13), the boundary flux for species  is then: where depends on the chosen units.For instance, in a 3D/2D system with volume concentration units of µM and surface concentration units of molecules/µm 2 , , where   is Avogadro's number.From Equation (S14), the surface reaction rate for  is: and the surface reaction rate for  is: In terms of these expressions, the system of PDEs is:

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  1 (Ω  ),  ∈ M of square-integrable functions on Ω  with square-integrable weak derivatives, and analogously for  1 (Γ  ),  ∈ , as well as the vector function spaces  1 (Ω  , R  ) =  1 (Ω)  .For solving the multi-domain reaction-diffusion equations, we introduce two product spaces  and  consisting of bulk fields and surface fields, respectively: To represent the solution fields ,  comprised of the separate bulk and surface components for the different species, we label

B.2 Discretization in time and space
SMART discretizes in time via an implicit (first-order) Euler scheme with time steps 0 =  0 <  1 < . . .<  −1 <   =  with timestep   =   − −1 .To discretize in space, we employ a conforming finite element mesh T defined relative to the parcellation of Ω into domains and surfaces, such that T  ⊆ T defines a submesh of Ω  for  ∈ M, and G  defines a (co-dimension one) mesh of Γ  for  ∈ Q.We define the finite element spaces of continuous piecewise linears P 1 defined relative to each (sub)mesh, and define the product spaces For each time step   , given discrete solutions  − ℎ and  − ℎ at the previous time  −1 , we then solve the nonlinear, coupled time-discrete problem: find for all  ∈  ℎ ,  ∈  ℎ , where    and    are defined by the forms  and  after implicit Euler time-discretization and depend on the given  − ℎ and  − ℎ .

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.

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,    defined over all mesh facets (faces for 2D parent mesh, edges for 3D parent mesh) and   defined over all mesh cells.For instance, in a 3D/2D model,   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   .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 , ,  or, in the case of surface species, the curvature .
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  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 , , , .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  (), flux values are updated as follows: 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 × 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"  is the multiplicative factor used to alter the time-step at the end of a given time step, i.e., the new time-step   is given by a in terms of the previous time-step   as   =    .
Table S1.Adaptive time stepping rules in SMART.

C.3 Axisymmetric models
SMART offers the ability to assume axisymmetry given a 2D mesh.Given a 2D mesh in the  −  plane, the symmetry axis is assumed to be  = 0.When this feature is turned on in SMART, the variational forms are adjusted accordingly, with the only necessary change being the multiplication by  in each integral (factors of 2 cancel throughout): 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.

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   and dephosphorylated through the cytosol with rate   41 .The SMART specifications are summarized in Tables S2 to S5.

Value/Equation Description
A_tot 1.000 µM Total cytosolic A k_kin 5.000 × 10 1 s −1 Rate constant for A phosphorylation VolSA 5.000 × 10 −2 µm Cytosolic volume to PM surface area ratio k_p 1.000 × 10 1 s −1 Rate constant for A dephosphorylation Given the two reactions summarized above, we can define the following flux terms in accordance with Appendix A.4: Accordingly the governing equation and boundary conditions are: 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  direction.A cell of thickness Δ with one region of membrane at  = 0 and the other at  = Δ is then predicted to have the following concentration profile at steady state: where

22/39
We ran these simulations for a thin slab where the thickness is ten times smaller than the length in the  and  directions, testing the effects of time-step refinement, mesh refinement, and changes in the diffusion coefficient.As an error metric, we computed the L 2 norm of the difference between the analytical steady state solution and our numerical solution at  = 1 s (  ℎ , ): Table S6.

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.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.The mesh statistics for each refinement are provided in Table S10.

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 where  and  are cylindrical coordinates.These implicit boundary representations were constructed by solving for  at each  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  value at the surface as a function of the cylindrical coordinate .Considering an implicitly (as above) or explicitly defined contour as a function of arc length  (), (), the scaled contour at a given value of  is  () (), ().We additionally require that this shape conserves cell volume, resulting in the constraint: For the circular contact region,  () = 1 and the above constraint is automatically satisfied.For a rectangular contact region with sides of length  and : resulting in the constraint  = .For our example, we chose  = √ 0.6 and  = √︁  0.6 for an aspect ratio of 0.6.For a star-shaped contact region: resulting in the constraint 2 2 0 +  2 1 = 2.For our example, we chose  0 = 0.98 and  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;  () for an ellipse with axes ,  is: and the associated constraint is  = 1.We choose the aspect ratio of the ellipse to be range(    ) range(    ) , where    and    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  functions that also depend on : where   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 µm.The nucleus shape remained unaltered across all simulations.

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.

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 µ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.

Figure 4 .Figure 4 .
Figure 4. Calcium dynamics within a realistic dendritic spine and calcium release unit.(cont.on next page)

Figure 5 .
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).

Figure 6 .
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 L 2 error between the computed and analytical solutions for all combinations of temporal and spatial refinement are shown in TableS6, TableS7and TableS8for diffusion coefficients  10, 100, and 1000 µm 2 s −1 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.

Figure 7 .
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.

) 18/ 39 Figure
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.

Figure S2 .
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.

Table S2 .
Compartments in phosphorylation model.Mesh statistics are given for the coarsest mesh (ℎ = 3.0)

Table S3 .
Species in phosphorylation model.

Table S5 .
Parameters in phosphorylation model.

Table S9 .
Refined meshes for numerical testing of the mechanotransduction model.

Table S10 .
Refined meshes for numerical testing of the dendritic spine model.

Table S15 .
Compartments in dendritic spine model.

Table S16 .
Species in dendritic spine model.

Table S19 .
Compartments in CRU model.

Table S20 .
Species in CRU model.
* * * 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.

Table S23 .
Compartments in ATP generation model.In simulations with cristae-localized membrane species, the IM was split into two regions as shown in Figure5B.The cristae portion had an area of 1.0197 µm 2 , with 28578 total surface elements.
* * Due to mass conservation,   is not treated as a separate variable, but solved for as  =   −  −  −   −  −   −   −   −

Table S26 .
Parameters for mitochondrial ATP generation.