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

Spatial modeling algorithms for reactions and transport (SMART) in biological cells | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Spatial modeling algorithms for reactions and transport (SMART) in biological cells View ORCID Profile Emmet A. Francis , View ORCID Profile Justin Laughlin , View ORCID Profile Jørgen S. Dokken , View ORCID Profile Henrik Finsberg , View ORCID Profile Christopher T. Lee , View ORCID Profile Marie E. Rognes , View ORCID Profile Padmini Rangamani doi: https://doi.org/10.1101/2024.05.23.595604 Emmet A. Francis 1 Department of Mechanical and Aerospace Engineering, University of California San Diego , La Jolla, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Emmet A. Francis Justin Laughlin 1 Department of Mechanical and Aerospace Engineering, University of California San Diego , La Jolla, CA, USA 2 Computational Engineering Division, Lawrence Livermore National Laboratory , Livermore, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Justin Laughlin Jørgen S. Dokken 3 Department of Numerical Analysis and Scientific Computing, Simula Research Laboratory , Oslo, Norway Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Jørgen S. Dokken Henrik Finsberg 4 Department of Computational Physiology, Simula Research Laboratory , Oslo, Norway Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Henrik Finsberg Christopher T. Lee 1 Department of Mechanical and Aerospace Engineering, University of California San Diego , La Jolla, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Christopher T. Lee Marie E. Rognes 3 Department of Numerical Analysis and Scientific Computing, Simula Research Laboratory , Oslo, Norway Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Marie E. Rognes Padmini Rangamani 1 Department of Mechanical and Aerospace Engineering, University of California San Diego , La Jolla, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Padmini Rangamani For correspondence: prangamani{at}ucsd.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF 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 – 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 – 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 – 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. Download figure Open in new tab 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 . Recent efforts have been made to accurately represent cellular and subcellular geometries with well-conditioned triangular and tetrahedral meshes 17 – 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 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. Download figure Open in new tab 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 (Ω m ) and 3 surface compartments (Γ 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 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 t final . Results are post-processed to examine changes in concentration over time and space. 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. 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. Download figure Open in new tab Figure 3. Spatial model of YAP/TAZ mechanotransduction for cells on micropatterned surfaces. 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. 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). 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 . Download figure Open in new tab Figure 4. Calcium dynamics within a realistic dendritic spine and calcium release unit. 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. 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 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. Download figure Open in new tab 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). 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 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 ℒ 2 norm of the error with respect to element size h ( Figure 6C ). Download figure Open in new tab 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 ℒ 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 µ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. 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 Δ 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 Table S11 ). Dendritic spine simulations were initially conducted on a mesh with h min = 0.009 µm and h max = 0.199 µm and a time step of Δ 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 µm and h max = 0.0383 µm, and our time step is adaptively adjusted from an initial value of Δ t = 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 106979 (18649), 855832 (146024), and 6846656 (1154871) 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 ( N , dimensions of the discrete solutions) are: 49194, 323328, 2299090; thus corresponding to a 6.6× and 7.1× 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 𝒪( N ) and 𝒪( 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 ). Download figure Open in new tab 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 – 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 – 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 (Ω) 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 Ω m ⊆ Ω are defined as a set of mesh cells with a common label or tag m ∈ ℳ, and similarly exterior boundaries and interior interfaces Γ q 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 Ω m and/or surface Γ 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 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 Ω m , and on or across the surfaces Γ 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 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 Ω m and the concentration of each species defined over each surface Γ 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 (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 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 + τ n was estimated as: 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: 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. 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 ℳ: index/label set for domains Ω m ( m ∈ ℳ). 𝒬: index/label set for the surfaces Γ q ( q ∈ 𝒬). ℐ, ℐ m , ℐ q : index/label sets for the species, the species in Ω m and species on Γ q , respectively. 𝒦, 𝒦 m , 𝒦 q , 𝒦 mq , 𝒦 mqn : index/label sets for the reactions, reactions in Ω m , reactions on Γ q , volume-surface reactions between species in Ω m and Γ q , and volume-surface-volume reactions between species in Ω m , Γ q , and Ω n . A.2 Multi-domain and multi-surface geometry representation We consider an open topologically D -dimensional manifold Ω ⊂ ℝ d for D ⩽ d = 1, 2, 3, and assume that Ω is partitioned into |ℳ| open and disjoint domains Ω m ⊂ ℝ d : each with (internal or external) boundary ∂ Ω m , and boundary ∂ Ω. In the SMART framework, each domain Ω m is referred to as a volume compartment. The boundary normal n m on ∂ Ω m is the outward pointing normal vector field to the boundary. The interface Γ mn between domains Ω m and Ω n is defined as the intersection of the closure of the domains: for m, n ∈ ℳ, and may or may not be empty. Each interface Γ mn may be further partitioned into one or more surfaces : where each Γ q is referred to as a surface compartment within SMART. We denote the total set of surfaces by 𝒬 = ∪ m,n ∈ℳ ∪ 𝒬 0 where 𝒬 0 also includes (sub)surfaces on the boundary of Ω. For each surface, we define the surface-to-neighbors map from the surface index q to its neighboring domain indices m, n : q ↦ m, n . If the surface is a real boundary surface, the map returns the single neighbor q ↦ m ,. Moreover, let time t ∈ [0, T ] for T > 0. Download figure Open in new tab 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. 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 Ω m by ℐ m for each (sub)domain Ω m . For each m , the species i ∈ ℐ m with concentrations for x ∈ Ω m and t ∈ (0, T ] satisfy reaction-transport equations of the form where ∂ t is the time-derivative, 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 . 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 i in domain Ω m which may be heterogeneous i.e., spatially varying and/or anisotropic i.e., tensor-valued. Remaining relative to Ω m , on (all) surfaces Γ q ⊆ Γ mn , we assume that the flux of species i is governed by a relation of the form: where defines surface fluxes, and the surface concentrations ν q are defined below. Surface species The surface concentrations , entering in (S6) above, are defined on Γ q ⊆ Γ mn , either as prescribed fields or via surface equations as follows: find ν q = ν q ( x, t ) for x ∈ Γ q , t ∈ (0, T ] such that where are surface reactions for each species j ∈ ℐ q . In the case of transport by surface diffusion, with surface diffusion coefficient for species j , surface gradient ∇ S and surface divergence ∇ S ·, we have: Boundary conditions For any (sub)boundary Γ q ⊂ ∂ Ω m such that Γ q ⊆ ∂ Ω, 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 ν q on any surface Γ q with non-empty boundary, zero flux boundary conditions are prescribed: where n q is the normal to the boundary line. Initial conditions Initial conditions are required for any for m ∈ ℳ, i ∈ ℐ m , and for q ∈ 𝒬, 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 Ω m , , the flux associated with a reactant or product species i for reaction k is given by: where 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 Γ q is , and the flux associated with a reactant or product species j for reaction k is given by: where 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 Ω m and Γ q , , the flux associated with a reactant or product volume species i for reaction k is given by: where 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 j that acts as a product or reactant, such reactions have an associated flux: where 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 Ω m , Γ q , and Ω n is 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 in Equation (S4) is given by summing the contribution from all volume reactions involving species i , : 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 i : 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 j , as well as any volume-surface or volume-surface-volume reactions involving species j : 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 Ω cyto that can bind to the surface species X on the plasma membrane domain Γ 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: From Equation (S13) , the boundary flux for species A 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 N A is Avogadro’s number. From Equation (S14) , the surface reaction rate for X is: and the surface reaction rate for B 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 H 1 (Ω m ), m ∈ ℳ of square-integrable functions on Ω m with square-integrable weak derivatives, and analogously for H 1 (Γ q ), q ∈ 𝒬, as well as the vector function spaces H 1 (Ω m , ℝ d ) = H 1 (Ω) 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: To represent the solution fields u, ν comprised of the separate bulk and surface components for the different species, we label 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 ∈ U and ν ∈ V such that for all ϕ ∈ U and ψ ∈ V : where both forms F and G are composed of sums over domains or surfaces and species: Furthermore, with the L 2 ( O )-inner product over any given domain O ⊂ Ω defined as we have defined and, for any q ∈ 𝒬 mn for given Ω m and Ω n interfacing Ω m via Γ q , finally: 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 < … < t n −1 < t n = T with timestep τ n = t n − t n −1 . To discretize in space, we employ a conforming finite element mesh 𝒯 defined relative to the parcellation of Ω into domains and surfaces, such that 𝒯 m ⊆ 𝒯 defines a submesh of Ω m for m ∈ ℳ, and 𝒢 q defines a (co-dimension one) mesh of Γ q for q ∈ 𝒢. We define the finite element spaces of continuous piecewise linears 𝒫 1 defined relative to each (sub)mesh, and define the product spaces For each time step t i , given discrete solutions and at the previous time t i −1 , we then solve the nonlinear, coupled time-discrete problem: find u h ∈ U h and ν h ∈ V h for all ϕ ∈ U h , ψ ∈ V h , where and are defined by the forms F and G 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, M f acet 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 curν . 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 curν 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: 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 τ new is given by a in terms of the previous time-step τ old as τ new = ζ τ old . View this table: View inline View popup Download powerpoint 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 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 π 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 k kin and dephosphorylated through the cytosol with rate k p 41 . The SMART specifications are summarized in Tables S2 to S5 . View this table: View inline View popup Download powerpoint Table S2. Compartments in phosphorylation model. Mesh statistics are given for the coarsest mesh ( h = 3.0) View this table: View inline View popup Download powerpoint Table S3. Species in phosphorylation model. View this table: View inline View popup Download powerpoint Table S4. Reactions in phosphorylation model. View this table: View inline View popup Download powerpoint Table S5. Parameters in phosphorylation model. 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 z direction. A cell of thickness Δ z with one region of membrane at z = 0 and the other at z = Δ z is then predicted to have the following concentration profile at steady state: 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 ℒ 2 norm of the difference between the analytical steady state solution and our numerical solution at : View this table: View inline View popup Download powerpoint Table S6. ℒ 2 error (in μM) between computed and analytical solution for different time steps and mesh sizes for D = 10.0 μm 2 s −1 . View this table: View inline View popup Download powerpoint Table S7. ℒ 2 error (in μM) between computed and analytical solution for different time steps and mesh sizes for D = 100.0 μm 2 s −1 . View this table: View inline View popup Download powerpoint Table S8. ℒ 2 error (in μM) between computed and analytical solution for different time steps and mesh sizes for D = 1000.0 μm 2 s −1 . 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 . View this table: View inline View popup Download powerpoint Table S9. Refined meshes for numerical testing of the mechanotransduction model. 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. Download figure Open in new tab 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 . View this table: View inline View popup Download powerpoint Table S10. Refined meshes for numerical testing of the dendritic spine model. 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 and the nuclear membrane surface satisfies 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 θ . 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 θ is r ( s ) T ( θ ), z ( s ). We additionally require that this shape conserves cell volume, resulting in the constraint: For the circular contact region, T ( θ ) = 1 and the above constraint is automatically satisfied. For a rectangular contact region with sides of length a and b : resulting in the constraint ab = π . For our example, we chose and for an aspect ratio of 0.6. For a star-shaped contact region: resulting in the constraint . 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 ( θ ) for an ellipse with axes c, d is: and the associated constraint is cd = 1. We choose the aspect ratio of the ellipse to be , 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 : 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 μ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. F.1 Mechanotransduction model View this table: View inline View popup Download powerpoint Table S11. Compartments in mechanotransduction model. View this table: View inline View popup Download powerpoint Table S12. Species in mechanotransduction model. View this table: View inline View popup Table S13. Reactions in mechanotransduction model. View this table: View inline View popup Table S14. Parameters in mechanotransduction model. Values were adopted directly from those reported in 33 F.2 Dendritic spine calcium model View this table: View inline View popup Download powerpoint Table S15. Compartments in dendritic spine model. View this table: View inline View popup Download powerpoint Table S16. Species in dendritic spine model. View this table: View inline View popup Table S17. Reactions in dendritic spine model. View this table: View inline View popup Download powerpoint Table S18. Parameters in dendritic spine model. F.3 Calcium release unit model View this table: View inline View popup Download powerpoint Table S19. Compartments in CRU model. View this table: View inline View popup Download powerpoint Table S20. Species in CRU model. View this table: View inline View popup Table S21. Reactions in CRU model. View this table: View inline View popup Table S22. Parameters in CRU model. F.4 ATP production in mitochondria View this table: View inline View popup Download powerpoint Table S23. Compartments in ATP generation model. View this table: View inline View popup Download powerpoint Table S24. Species in ATP generation model. View this table: View inline View popup Download powerpoint Table S25. Reactions in ATP generation model. View this table: View inline View popup Table S26. Parameters for mitochondrial ATP generation. 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. 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. Footnotes ↵ * meg{at}simula.no , prangamani{at}ucsd.edu https://zenodo.org/doi/10.5281/zenodo.10019463 https://zenodo.org/doi/10.5281/zenodo.10480303 https://zenodo.org/doi/10.5281/zenodo.11268944 https://zenodo.org/doi/10.5281/zenodo.11252054 References 1. ↵ Yue , R. & Dutta , A. Computational systems biology in disease modeling and control, review and perspectives . npj Syst. Biol. Appl . 8 , 1 – 16 , DOI: 10.1038/s41540-022-00247-4 ( 2022 ). OpenUrl CrossRef 2. ↵ Mogilner , A. , Wollman , R. & Marshall , W. F. Quantitative Modeling in Cell Biology: What Is It Good for? Dev . Cell 11 , 279 – 287 , DOI: 10.1016/j.devcel.2006.08.004 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 3. ↵ Kitano , H. Computational systems biology . Nature 420 , 206 – 210 , DOI: 10.1038/nature01254 ( 2002 ). OpenUrl CrossRef PubMed Web of Science 4. ↵ Wainwright , S. A. Form and Function in Organisms . Am. Zool . 28 , 671 – 680 , DOI: 10.1093/icb/28.2.671 ( 1988 ). OpenUrl CrossRef 5. ↵ Herman , M. A. et al. A Unifying Framework for Understanding Biological Structures and Functions Across Levels of Biological Organization . Integr. Comp. Biol . 61 , 2038 – 2047 , DOI: 10.1093/icb/icab167 ( 2022 ). OpenUrl CrossRef 6. ↵ Peddie , C. J. et al. Volume electron microscopy . Nat. Rev. Methods Primers 2 , 1 – 23 , DOI: 10.1038/s43586-022-00131-9 ( 2022 ). OpenUrl CrossRef 7. Villinger , C. et al. FIB/SEM tomography with TEM-like resolution for 3D imaging of high-pressure frozen cells . Histochem. Cell Biol . 138 , 549 – 556 , DOI: 10.1007/s00418-012-1020-6 ( 2012 ). OpenUrl CrossRef PubMed 8. ↵ Heinrich , L. et al. Whole-cell organelle segmentation in volume electron microscopy . Nature 599 , 141 – 146 , DOI: 10.1038/s41586-021-03977-3 ( 2021 ). OpenUrl CrossRef 9. ↵ McCafferty , C. L. et al. Integrating cellular electron microscopy with multimodal data to explore biology across space and time . Cell 187 , 563 – 584 , DOI: 10.1016/j.cell.2024.01.005 ( 2024 ). OpenUrl CrossRef 10. ↵ Stone , M. B. , Shelby , S. A. & Veatch , S. L. Super-Resolution Microscopy: Shedding Light on the Cellular Plasma Membrane . Chem. Rev . 117 , 7457 – 7477 , DOI: 10.1021/acs.chemrev.6b00716 ( 2017 ). OpenUrl CrossRef PubMed 11. ↵ Schermelleh , L. et al. Super-resolution microscopy demystified . Nat. Cell Biol . 21 , 72 – 84 , DOI: 10.1038/s41556-018-0251-8 ( 2019 ). OpenUrl CrossRef PubMed 12. ↵ Hake , J. et al. Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit: Calcium sparks in reconstructed release unit . The J. Physiol . 590 , 4403 – 4422 , DOI: 10.1113/jphysiol.2012.227926 ( 2012 ). OpenUrl CrossRef PubMed 13. Bell , M. K. , Holst , M. V. , Lee , C. T. & Rangamani , P. Dendritic spine morphology regulates calcium-dependent synaptic weight change . The J. Gen. Physiol . 154 , e202112980 , DOI: 10.1085/jgp.202112980 ( 2022 ). OpenUrl CrossRef 14. ↵ Garcia , G. C. et al. Mitochondrial morphology provides a mechanism for energy buffering at synapses . Sci. Reports 9 , 18306 , DOI: 10.1038/s41598-019-54159-1 ( 2019 ). OpenUrl CrossRef 15. ↵ Holash , R. J. & MacIntosh , B. R. A stochastic simulation of skeletal muscle calcium transients in a structurally realistic sarcomere model using MCell . PLoS Comput. Biol . 15 , e1006712 , DOI: 10.1371/journal.pcbi.1006712 ( 2019 ). OpenUrl CrossRef 16. ↵ Daversin-Catty , C. , Richardson , C. N. , Ellingsrud , A. J. & Rognes , M. E. Abstractions and Automated Algorithms for Mixed Domain Finite Element Methods . ACM Transactions on Math. Softw . 47 , 31 :1–31:36 , DOI: 10.1145/3471138 ( 2021 ). OpenUrl CrossRef 17. ↵ Wang , X. & Danuser , G. Remeshing flexible membranes under the control of free energy . PLOS Comput. Biol . 18 , e1010766 , DOI: 10.1371/journal.pcbi.1010766 ( 2022 ). OpenUrl CrossRef 18. Means , S. et al. Reaction Diffusion Modeling of Calcium Dynamics with Realistic ER Geometry . Biophys. J . 91 , 537 – 557 , DOI: 10.1529/biophysj.105.075036 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 19. ↵ Lee , C. T. et al. 3D mesh processing using GAMer 2 to enable reaction-diffusion simulations in realistic cellular geometries . PLOS Comput. Biol . 16 , e1007756 , DOI: 10.1371/journal.pcbi.1007756 ( 2020 ). OpenUrl CrossRef 20. ↵ Lee , C. T. et al. An Open-Source Mesh Generation Platform for Biophysical Modeling Using Realistic Cellular Geometries . Biophys. J . 118 , 1003 – 1008 , DOI: 10.1016/j.bpj.2019.11.3400 ( 2020 ). OpenUrl CrossRef 21. ↵ Venkatraman , K. et al. Cristae formation is a mechanical buckling event controlled by the inner mitochondrial membrane lipidome . The EMBO J . 42 , e114054 , DOI: 10.15252/embj.2023114054 ( 2023 ). OpenUrl CrossRef 22. ↵ Mesa , M. H. , Garcia , G. C. , Hoerndli , F. J. , McCabe , K. J. & Rangamani , P. Spine apparatus modulates Ca2+ in spines through spatial localization of sources and sinks , DOI: 10.1101/2023.09.22.558941 ( 2023 ). OpenUrl Abstract / FREE Full Text 23. ↵ Alnæs , M. et al. The FEniCS Project Version 1.5 . Archive Numer. Softw . Vol 3 , Starting Point and Frequency: Year: 2013 , DOI: 10.11588/ANS.2015.100.20553 ( 2015 ). OpenUrl CrossRef 24. ↵ Laughlin , J. G. et al. SMART: Spatial Modeling Algorithms for Reactions and Transport . J. Open Source Softw . 8 , 5580 , DOI: 10.21105/joss.05580 ( 2023 ). OpenUrl CrossRef 25. ↵ Laughlin , J. et al. RangamaniLabUCSD/smart: V2.2.1. Zenodo , DOI: 10.5281/zenodo.11268687 ( 2024 ). OpenUrl CrossRef 26. ↵ Schaff , J. , Fink , C. , Slepchenko , B. , Carson , J. & Loew , L. A general computational framework for modeling cellular structure and function . Biophys. J . 73 , 1135 – 1146 , DOI: 10.1016/S0006-3495(97)78146-3 ( 1997 ). OpenUrl CrossRef PubMed Web of Science 27. ↵ Cowan , A. E. , Moraru , I. I. , Schaff , J. C. , Slepchenko , B. M. & Loew , L. M. Spatial Modeling of Cell Signaling Networks . In Methods in Cell Biology , vol. 110 , 195 – 221 , DOI: 10.1016/B978-0-12-388403-9.00008-4 ( Elsevier , 2012 ). OpenUrl CrossRef PubMed Web of Science 28. ↵ Kerr , R. A. et al. Fast Monte Carlo Simulation Methods for Biological Reaction-Diffusion Systems in Solution and on Surfaces . SIAM J. on Sci. Comput . 30 , 3126 – 3149 , DOI: 10.1137/070692017 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 29. ↵ Hucka , M. et al. The systems biology markup language (SBML): A medium forrepresentation and exchange of biochemical network models . Bioinformatics 19 , 524 – 531 , DOI: 10.1093/bioinformatics/btg015 ( 2003 ). OpenUrl CrossRef PubMed Web of Science 30. ↵ Wu , Y. et al. Contacts between the endoplasmic reticulum and other membranes in neurons . Proc. Natl. Acad. Sci. United States Am . 114 , E4859 – E4867 , DOI: 10.1073/pnas.1701078114 ( 2017 ). OpenUrl Abstract / FREE Full Text 31. ↵ Balay , S. , Gropp , W. , McInnes , L. C. & Smith , B. F. PETSc, the portable, extensible toolkit for scientific computation ( 1998 ). 32. ↵ Sun , M. , Spill , F. & Zaman , M. H. A Computational Model of YAP/TAZ Mechanosensing . Biophys. J . 110 , 2540 – 2550 , DOI: 10.1016/j.bpj.2016.04.040 ( 2016 ). OpenUrl CrossRef 33. ↵ Scott , K. E. , Fraley , S. I. & Rangamani , P. A spatial model of YAP/TAZ signaling reveals how stiffness, dimensionality, and shape contribute to emergent outcomes . Proc. Natl. Acad. Sci . 118 , e2021571118 , DOI: 10.1073/pnas.2021571118 ( 2021 ). OpenUrl Abstract / FREE Full Text 34. ↵ Maduram , J. H. , Goluch , E. , Hu , H. , Liu , C. & Mrksich , M. Subcellular curvature at the perimeter of micropatterned cells influences lamellipodial distribution and cell polarity . Cell motility cytoskeleton 65 , 841 – 852 , DOI: 10.1002/cm.20305 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 35. ↵ Bell , M. , Bartol , T. , Sejnowski , T. & Rangamani , P. Dendritic spine geometry and spine apparatus organization govern the spatiotemporal dynamics of calcium . J. Gen. Physiol . 151 , 1017 – 1034 , DOI: 10.1085/jgp.201812261 ( 2019 ). OpenUrl Abstract / FREE Full Text 36. ↵ Leung , A. , Ohadi , D. , Pekkurnaz , G. & Rangamani , P. Systems modeling predicts that mitochondria ER contact sites regulate the postsynaptic energy landscape . NPJ Syst. Biol. Appl . 7 , 26 , DOI: 10.1038/s41540-021-00185-7 ( 2021 ). OpenUrl CrossRef 37. ↵ Parinandi , N. L. & Hund , T. J. Li , J. , Richmond , B. & Hong , T. Organization of Ca2+ Signaling Microdomains in Cardiac Myocytes . In Parinandi , N. L. & Hund , T. J. (eds.) Cardiovascular Signaling in Health and Disease ( Springer , Cham (CH ), 2022 ). 38. ↵ Hoshijima , M. et al. CCDB:3603, MUS MUSCULUS, T-tubules, sarcoplasmic reticulum, myocyte , DOI : doi: 10.7295/W9CCDB3603 ( 2004 ). OpenUrl CrossRef 39. ↵ Mendelsohn , R. et al. Morphological principles of neuronal mitochondria . The J. Comp. Neurol . 530 , 886 – 902 , DOI: 10.1002/cne.25254 ( 2022 ). OpenUrl CrossRef 40. ↵ Garcia , G. C. , Gupta , K. , Bartol , T. M. , Sejnowski , T. J. & Rangamani , P. Mitochondrial morphology governs ATP production rate . J. Gen. Physiol . 155 , DOI: 10.1085/jgp.202213263 ( 2023 ). OpenUrl CrossRef 41. ↵ Meyers , J. , Craig , J. & Odde , D. J. Potential for Control of Signaling Pathways via Cell Size and Shape . Curr. Biol . 16 , 1685 – 1693 , DOI: 10.1016/j.cub.2006.07.056 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 42. ↵ Rangamani , P. et al. Decoding Information in Cell Shape . Cell 154 , 1356 – 1369 , DOI: 10.1016/j.cell.2013.08.026 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 43. Neves , S. R. et al. Cell Shape and Negative Links in Regulatory Motifs Together Control Spatial Information Flow in Signaling Networks . Cell 133 , 666 – 680 , DOI: 10.1016/j.cell.2008.04.025 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 44. ↵ Calizo , R. C. et al. Cell shape regulates subcellular organelle location to control early Ca2+ signal dynamics in vascular smooth muscle cells . Sci. Reports 10 , 17866 , DOI: 10.1038/s41598-020-74700-x ( 2020 ). OpenUrl CrossRef 45. ↵ Johnson , G. T. et al. Building the next generation of virtual cells to understand cellular biology . Biophys. J . 122 , 3560 – 3569 , DOI: 10.1016/j.bpj.2023.04.006 ( 2023 ). OpenUrl CrossRef 46. ↵ Viana , M. P. et al. Integrated intracellular organization and its variations in human iPS cells . Nature 613 , 345 – 354 , DOI: 10.1038/s41586-022-05563-7 ( 2023 ). OpenUrl CrossRef 47. ↵ Xu , C. S. et al. An open-access volume electron microscopy atlas of whole cells and tissues . Nature 599 , 147 – 151 , DOI: 10.1038/s41586-021-03992-4 ( 2021 ). OpenUrl CrossRef 48. ↵ Edwards , J. et al. VolRoverN: Enhancing Surface and Volumetric Reconstruction for Realistic Dynamical Simulation of Cellular and Subcellular Function . Neuroinformatics 12 , 277 – 289 , DOI: 10.1007/s12021-013-9205-2 ( 2014 ). OpenUrl CrossRef PubMed 49. ↵ Hu , Y. , Schneider , T. , Wang , B. , Zorin , D. & Panozzo , D. Fast tetrahedral meshing in the wild . ACM Transactions on Graph . 39 , 117 :117:1–117:117:18 , DOI: 10.1145/3386569.3392385 ( 2020 ). OpenUrl CrossRef 50. ↵ Mitusch , S. K. , Funke , S. W. & Dokken , J. S. Dolfin-adjoint 2018.1: Automated adjoints for FEniCS and Firedrake . J. Open Source Softw . 4 , 1292 , DOI: 10.21105/joss.01292 ( 2019 ). OpenUrl CrossRef 51. ↵ Novak , I. L. et al. Diffusion on a Curved Surface Coupled to Diffusion in the Volume: Application to Cell Biology . J. computational physics 226 , 1271 – 1290 , DOI: 10.1016/j.jcp.2007.05.025 ( 2007 ). OpenUrl CrossRef PubMed 52. ↵ Hepburn , I. , Chen , W. , Wils , S. & De Schutter , E. STEPS: Efficient simulation of stochastic reaction–diffusion models in realistic morphologies . BMC Syst. Biol . 6 , 36 , DOI: 10.1186/1752-0509-6-36 ( 2012 ). OpenUrl CrossRef PubMed 53. ↵ Hepburn , I. et al. Vesicle and reaction-diffusion hybrid modeling with STEPS . Commun. Biol . 7 , 1 – 15 , DOI: 10.1038/s42003-024-06276-5 ( 2024 ). OpenUrl CrossRef 54. ↵ Roberts , E. , Stone , J. E. & Luthey-Schulten , Z. Lattice microbes: High-performance stochastic simulation method for the reaction-diffusion master equation . J. Comput. Chem . 34 , 245 – 255 , DOI: 10.1002/jcc.23130 ( 2013 ). OpenUrl CrossRef PubMed 55. ↵ Hoops , S. et al. COPASI—a COmplex PAthway SImulator . Bioinformatics 22 , 3067 – 3074 , DOI: 10.1093/bioinformatics/btl485 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 56. ↵ Clark , A. P. , Chowkwale , M. , Paap , A. , Dang , S. & Saucerman , J. J. Logic-based modeling of biological networks with Netflux . bioRxiv 2024.01.11.575227 , DOI: 10.1101/2024.01.11.575227 ( 2024 ). OpenUrl Abstract / FREE Full Text 57. ↵ Geuzaine , C. & Remacle , J.-F. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities: THE GMSH PAPER . Int. J. for Numer. Methods Eng . 79 , 1309 – 1331 , DOI: 10.1002/nme.2579 ( 2009 ). OpenUrl CrossRef 58. ↵ Francis , E. A. , Finsberg , H. N. & Dokken , J. S. Biological test cases implemented in SMART . Zenodo , DOI: 10.5281/zenodo.11268945 ( 2024 ). OpenUrl CrossRef 59. ↵ Alnæs , M. S. , Logg , A. , Ølgaard , K. B. , Rognes , M. E. & Wells , G. N. Unified form language: A domain-specific language for weak formulations of partial differential equations . ACM Transactions on Math. Softw . 40 , 9 :1–9:37 , DOI: 10.1145/2566630 ( 2014 ). OpenUrl CrossRef 60. ↵ Grecco , H. E. Pint: Physical quantities module . 61. ↵ Dalcin , L. D. , Paz , R. R. , Kler , P. A. & Cosimo , A. Parallel distributed computing using Python . Adv. Water Resour . 34 , 1124 – 1139 , DOI: 10.1016/j.advwatres.2011.04.013 ( 2011 ). OpenUrl CrossRef GeoRef 62. ↵ Meurer , A. et al. SymPy: Symbolic computing in Python . PeerJ Comput. Sci . 3 , e103 , DOI: 10.7717/peerj-cs.103 ( 2017 ). OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted May 25, 2024. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Spatial modeling algorithms for reactions and transport (SMART) in biological cells Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Spatial modeling algorithms for reactions and transport (SMART) in biological cells Emmet A. Francis , Justin Laughlin , Jørgen S. Dokken , Henrik Finsberg , Christopher T. Lee , Marie E. Rognes , Padmini Rangamani bioRxiv 2024.05.23.595604; doi: https://doi.org/10.1101/2024.05.23.595604 Share This Article: Copy Citation Tools Spatial modeling algorithms for reactions and transport (SMART) in biological cells Emmet A. Francis , Justin Laughlin , Jørgen S. Dokken , Henrik Finsberg , Christopher T. Lee , Marie E. Rognes , Padmini Rangamani bioRxiv 2024.05.23.595604; doi: https://doi.org/10.1101/2024.05.23.595604 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Systems Biology Subject Areas All Articles Animal Behavior and Cognition (7651) Biochemistry (17746) Bioengineering (13928) Bioinformatics (42066) Biophysics (21499) Cancer Biology (18650) Cell Biology (25579) Clinical Trials (138) Developmental Biology (13409) Ecology (19947) Epidemiology (2067) Evolutionary Biology (24374) Genetics (15633) Genomics (22557) Immunology (17775) Microbiology (40505) Molecular Biology (17217) Neuroscience (88796) Paleontology (667) Pathology (2845) Pharmacology and Toxicology (4836) Physiology (7664) Plant Biology (15179) Scientific Communication and Education (2047) Synthetic Biology (4304) Systems Biology (9839) Zoology (2272)