Finite element analysis predicts Ca2+ microdomains within tubular-sarcoplasmic reticular junctions of amphibian skeletal muscle

A finite element analysis modelled diffusional generation of steady-state Ca2+ microdomains within skeletal muscle transverse (T)-tubular-sarcoplasmic reticular (SR) junctions, sites of ryanodine receptor (RyR)-mediated SR Ca2+ release. It used established quantifications of sarcomere and T-SR anatomy (radial diameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=220 \, \mathrm{n}\mathrm{m}$$\end{document}d=220nm; axial distance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w=12 \, \mathrm{n}\mathrm{m}$$\end{document}w=12nm). Its boundary SR Ca2+ influx densities,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J}_{\mathrm{influx}}$$\end{document}Jinflux, reflected step impositions of influxes, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\it {\Phi }_{\mathrm{influx}}={J}_{\mathrm{influx}}\left(\frac{\pi {d}^{2}}{4}\right),$$\end{document}Φinflux=Jinfluxπd24, deduced from previously measured Ca2+ signals following muscle fibre depolarization. Predicted steady-state T-SR junctional edge [Ca2+], [Ca2+]edge, matched reported corresponding experimental cytosolic [Ca2+] elevations given diffusional boundary efflux \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\it \it {\Phi }_{\mathrm{efflux}}=\frac{D [ {{{\mathrm{Ca}}^{2+}}}]_{\mathrm{edge}}}{\lambda } (\pi dw),$$\end{document}Φefflux=D[Ca2+]edgeλ(πdw), established cytosolic Ca2+ diffusion coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(D = 4 \times {10}^{7} \mathrm{nm}^{2}/\mathrm{s})$$\end{document}(D=4×107nm2/s) and exit length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda = 9.2 \, \mathrm{n}\mathrm{m}$$\end{document}λ=9.2nm. Dependences of predicted [Ca2+]edge upon \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J}_{\mathrm{influx}}$$\end{document}Jinflux then matched those of experimental [Ca2+] upon Ca2+ release through their entire test voltage range. The resulting model consistently predicted elevated steady-state T-SR junctional ~ µM-[Ca2+] elevations radially declining from maxima at the T-SR junction centre along the entire axial T-SR distance. These [Ca2+] heterogeneities persisted through 104- and fivefold, variations in D and w around, and fivefold reductions in d below, control values, and through reported resting muscle cytosolic [Ca2+] values, whilst preserving the flux conservation (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\it \it {\Phi }_{\mathrm{influx}}={\Phi }_{\mathrm{efflux}})$$\end{document}Φinflux=Φefflux) condition, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\left[\mathrm{C}{\mathrm{a}}^{2+}\right]}_{\mathrm{edge}}=\frac{\lambda {dJ}_{\mathrm{influx}}}{4Dw}$$\end{document}Ca2+edge=λdJinflux4Dw. Skeletal muscle thus potentially forms physiologically significant ~ µM-[Ca2+] T-SR microdomains that could regulate cytosolic and membrane signalling molecules including calmodulin and RyR, These findings directly fulfil recent experimental predictions invoking such Ca2+ microdomains in observed regulatory effects upon Na+ channel function, in a mechanism potentially occurring in similar restricted intracellular spaces in other cell types.

] T-SR microdomains that could regulate cytosolic and membrane signalling molecules including calmodulin and RyR, These findings directly fulfil recent experimental predictions invoking such Ca 2+ microdomains in observed regulatory effects upon Na + channel function, in a mechanism potentially occurring in similar restricted intracellular spaces in other cell types.
Intracellular endoplasmic or sarcoplasmic reticular (SR) membrane systems gating store Ca 2+ release into the cytosol following surface membrane activation, often involving ryanodine receptor (RyR) activation, occur widely amongst cell types. These include both excitable (skeletal, cardiac and smooth muscle, and cerebellar Purkinje 1-3 , hippocampal 4 and other central nervous system neurones 5,6 ) and non-excitable, including thrombocyte, cell types 7,8 . These intracellular membranes often form appositions with surface membrane with proximities (< 10-30 nm) permitting direct protein-protein/lipid interaction 9,10 though not accommodating entire organelles. Their intervening electron-dense cytosol could also reflect local concentrations of proteins, lipids or ions. Electron microscopic sections can reveal parallel alignments extending over ~ 100-400 nm distances without fusion of the component membranes potentially offering restricted diffusion spaces permitting ion, including Ca 2+ , accumulation and microdomain formation.
In skeletal and cardiac muscle, following surface membrane propagation, Na + channel mediated action potentials are conducted into the cellular interior at regular intervals along the muscle length through electrically continuous transverse (T-) tubular membranes. At specific regions, these come geometrically close (~ 12 nm) to, whilst remaining electrically isolated from, terminal cisternal membranes of the SR Ca 2+ store. The resulting T-SR triad and dyad junctions are strategic to excitation-contraction coupling [11][12][13] . In cardiac muscle, tubular
These paradoxical findings prompted suggestions that RyR1-mediated Ca 2+ release took place into a microdomain in the vicinity of both the SR RyR1 and the T-tubular membrane Na v 1.4 and that the local elevation in Ca 2+ concentration, [Ca 2+ ] TSR , would then modify Na v 1.4 function 30 . Such a hypothesis would predict contrasting increases and decreases in local microdomain [Ca 2+ ] following caffeine, and dantrolene or CPA challenge. The narrow, ~ 12 nm T-SR junctions that could form spaces with restricted intracellular diffusion close to the RyR1 Ca 2+ -release sites might be implicated in such microdomain formation. These could result in changes in local in vivo [Ca 2+ ] TSR , distinct from those of [Ca 2+ ] i in the remaining bulk cytosol. This could explain the contrasting actions of RyR1 agonists and RyR1 or SERCA antagonists on I Na through correspondingly contrasting effects on local Ca 2+ or Ca 2+ /CaM levels, to which the Na v 1.4 would be directly or indirectly exposed, despite their similar effects on bulk cytosolic [Ca 2+ ] i 33,34 . Direct experimental explorations for such [Ca 2+ ] TSR microdomains possibly using fluorescent Ca 2+ indicator methods need to address the small dimensions and dispersed nature of the T-SR compartment [35][36][37] . The present complementary approach applies diffusional modelling techniques 38 to explore the physical parameters permitting accumulation or depletion of released SR Ca 2+ within the T-SR junction. It demonstrated that established anatomical and physiological features related to skeletal muscle excitation-contraction coupling are physically compatible with generation of significant Ca 2+ microdomains in both activated and resting muscle fibres. We then discuss their possible physiological effects both in myocyte T-SR junctions and in similar membrane appositions in other cell types.

Results
T-SR junction structure represented using a formal geometric model. We employed anatomical, optical and electron microscope quantifications of sarcomere and T-SR junction structure from amphibian twitch fibres as the muscle type for which the fullest data are available [39][40][41][42][43][44] . This provided the required details of T-tubular-sarcoplasmic reticular (T-SR) junction anatomy for the modelling studies (Table 1). First, the reported values of sarcomere length l, fibre diameter a, and relative tubular (T) to surface (S) membrane area reflected in the ratio of their respective capacitances C T /C S , yielded the sarcomeric surface membrane area: and the total tubular membrane surface area: Morphometric electron micrographic estimates of the proportion ξ of the total tubular membrane area A T accounted for by T-SR junctions 44 then gave the total T-tubular membrane area contributing to T-SR junction structures A TSR = ξA T .,w,d of SR membrane enclosed in a single T-SR junction. The latter was accordingly modelled as a circularly symmetrical disk-shaped space of dimensions w = 12 nm and d = 220 nm (Fig. 1, Ca 2+ diffusion into and through a single T-SR junction modelled by finite element analysis. The finite element analysis of steady state Ca 2+ diffusion through a single T-SR space used the circular ends of the geometry defined above to represent its respective T-tubular and SR membranes. The rim separating their edges connected the T-SR and remaining intracellular spaces. Influx boundary conditions were supplied by a steadystate and uniform Ca 2+ influx density J influx across the SR membrane face of each individual T-SR junction: In contrast, the T-tubular face represented a zero-flux boundary surface. The Ca 2+ then diffuses through and leaves the T-SR space at the rim with diffusion coefficient D: Away from the boundaries where p = 0: This gives efflux equation: Its Fick's Law constant D/λ comprises the Ca 2+ diffusion coefficient D and empirical exit length λ. The latter provides a geometrical parametrization of the re-uptake of the dissipated Ca 2+ from a well stirred cytosolic compartment by SERCA activity without saturation. The λ term represents the only free parameter in the entire modelling analysis. Diffusional processes within the T-SR junction were represented by superimposing a finite element mesh upon the T-SR junction geometry using PDE Toolbox, dividing that 3D geometry into fine tetrahedral elements of specified maximum length χ (Fig. 1, Right panels A-C). Different mesh sizes progressively divided the T-SR space into (A) 12 nm, (B) 6 nm and (C) 3 nm tetrahedral elements; finer mesh sizes were used where modelling investigated axial in addition to radial [Ca 2+ ] gradients and in duplicate runs matching different mesh sizes to validate those used in the reconstructions. The FEM and MATLAB solved the equations specified for the system and specified input parameters producing a [Ca 2+ ] dataset in the form of an array, whose spatial resolution was determined by the fineness of the mesh, set by the maximum element length χ. Its steady state solutions satisfied the overall conservation condition, Figure 1. T-SR junction structure represented using a formal geometric model divided into finite elements. Left: Geometrical representation of formalized T-SR junction represented as two disks, respectively of tubular (T) and sarcoplasmic reticular (SR) membranes, of diameter, d = 220 nm, within the radial (xy) plane, separated along an axial (z) direction by a T-SR distance, w = 12 nm. Right: Superimposition of finite elements dividing T-SR junction geometry into tetrahedral 'elements' with specified maximum lengths of (A) 12 nm, or 100%, (B) 6 nm or 50%, and (C) 3 nm or 25%, of the T-SR distance respectively. Reduced maximum tetrahedral lengths increase the number of elements, within which each time-dependent solution is generated, increasing spatial resolution in the solution concentration profile. This condition was used in checks for the steady state condition in detailed explorations of the effect of specific parameters that follow.
Sarcoplasmic reticulum Ca 2+ release producing Ca 2+ microdomains characterized by radial concentration gradients in the T-SR junction. The modelling process first sought three-dimensional (3D) reconstructions of radial steady state [Ca 2+ ] distributions resulting from diffusional processes using the computational T-SR parameters listed in Table 1. Ca 2+ release influx densities into each individual T-SR junction, J influx, used in the MATLAB program PDE Toolbox, could be related to previously reported 45  Antipyrylazo III absorbances in amphibian skeletal muscle fibres subject to voltage clamp steps from a − 90 mV resting to a 0 mV test membrane potential 45 reported a value of d[Ca 2+ ]/dt = 180 μmol/(dm 3 s). The Ca 2+ then diffuses through the T-SR space down its resulting concentration gradients with the diffusion constant D = 4.0 × 10 7 nm 2 /s previously reported for amphibian skeletal muscle [46][47][48] . The Ca 2+ finally leaves the T-SR junction space effluxing into surrounding cytosol at the edge of the T-SR junction across diffusional area πdw at a rate driven by [Ca 2+ ] TSR edge . This proved close and proportional to reported experimental peak cytosolic Ca 2+ concentration, [Ca 2+ ] max at the explored 0 mV test voltage with the use of an exit length value λ = 9.2 nm 45 .
An overall rate constant describing the dependence of the summed Ca 2+ fluxes upon [Ca 2+ ] TSR edge could be determined using the predicted number of T-SR junctions in unit muscle volume, The total Ca 2+ efflux into unit muscle volume is then: The constant of proportionality describing this linear dependence on [Ca 2+ ] TSR edge is then.
This resulting rate constant is smaller than but comparable to experimental rate constants describing eventual SR resequestration of the released cytosolic Ca 2+ . Thus, previous experimental studies suggested rate constants for such unsaturable SR Ca 2+-ATPase mediated Ca 2+ uptake around 22.3 ± 8.14/s under similarly steady state conditions where Ca 2+ binding to remaining, saturable, fast-exchanging cytosolic binding sites was constant 49 .  Figure 2E demonstrates an axially sliced Ca 2+ microdomain extending to the edge of the T-SR junction with [Ca 2+ ] radially falling ~ fivefold over ~ 100 nm from the T-SR junction centre. Figure 2F   www.nature.com/scientificreports/ was assigned a number referenced using a command retrieving the solution from the node closest to a specified 3-dimensional Cartesian co-ordinate (https:// www. mathw orks. com/ help/ pde/ ug/ heat-trans fer-probl em-withtempe rature-depen dent-prope rties. html). The solution from the node closest to selected points representing the edge and centre of the geometry, and half-way between these was plotted against time through the time trajectory of the computation (Fig. 2G). These typically reached steady state over a ~ 0.2 ms exponential timecourse following the onset of the imposed J influx, a timescale ~ 1-2 orders of magnitude shorter than experimentally measured Ca 2+ transients in vivo 45 . The resulting steady state [Ca 2+ ] at axial distances within the T-SR junction could additionally be obtained by interpolation for plotting against radial distance from the centre of the T-SR junction (Fig. 2H). Comparing [Ca 2+ ] timecourses with imposition and termination of the Ca 2+ fluxes confirmed that [Ca 2+ ] recovered back to its initial starting value, as expected for a first-order diffusional system (Fig. 2I). Finally, solutions obtained using varying 6 nm and 3 nm mesh sizes gave [Ca 2+ ] microdomain characteristics in close agreement, validating the computational parameters used in our finite element analysis (Supplementary  Table S1).

Persistent Ca 2+ microdomains with [Ca 2+ ] graded with Ca 2+ flux densities through varied test voltages.
The T-SR junction model was then tested against varying J influx holding constant the remaining diffusional, D and λ, and T-SR geometrical parameters, d and w (Supplementary Table S2 45 . Firstly, all the test voltages were associated with T-SR junction Ca 2+ microdomains (Fig. 3). Figure 3A plots diffusion equation solutions with time after imposing the Ca 2+ influx to reach steady state [Ca 2+ ] values at the centre and edge of the T-SR junction, and halfway between these points. Figure 3B shows these dependences of the steady state [Ca 2+ ] with radial distance close to the T-tubular (superscript 'T') and SR membranes (superscript 'SR'), and halfway along their axial distance (superscript 'TSR'). [Ca 2+ ] declined from its maximum values with radial distance from the T-SR junction centre (subscripts below: 'centre') towards the T-SR edge (subscript 'edge') with similar readouts close to the T-tubular (subscript: 'T') and SR membranes (subscript: 'SR'), and half-way (subscript: '50%') along their axial distance.
Secondly, quantifications of Ca 2+ microdomains arising from the range of J influx values and corresponding test voltages showed that these were similar in form, consistent with the linearity expected of diffusional processes. This emerged from comparing ratios between [Ca 2+ ] at the T-SR centre with those at the edge and halfway between (50%) at the T-tubular and SR membranes, and within the T-SR space.  45 . Figure 3D shows a linear relationship between [Ca 2+ ] edge , and   Ca 2+ microdomains with varied Ca 2+ diffusion coefficients. Skeletal muscle T-SR junction properties vary with physiological conditions, both within and between muscle fibres, individual muscles and muscle types. Furthermore, such surface-cytoplasmic membrane appositions also occur in and vary amongst other cell types. Nevertheless, the Ca 2+ microdomains robustly persisted with wide variations involving previously experimentally reported Ca 2+ diffusion constants, D, and T-SR junction geometries represented by their axial distances w, and diameters, d. These further computations varied each parameter in turn holding the remaining variables constant under conditions of fixed J influx .
Firstly, we extended the modelling beyond the reported skeletal muscle cytosolic Ca 2+ diffusion coefficient D = 4.0 × 10 7 nm 2 /s [46][47][48] . Higher reported values reach 5.2 × 10 8 nm 2 /s in other cell types 50 (cf: [51][52][53][54] ) and 1 × 10 9 nm 2 s at infinite dilution in vitro 55,56 (Supplementary Table S4). Lower D values might result from local concentrations of proteins, lipids or ions, suggested by the reported electron-dense T-SR junction cytosol 44 . Solutions with D between 4.0 × 10 5 and 1 × 10 9 nm 2 /s at constant w, d and J influx , continued to converge with correspondingly varied relative timescales (Fig. 4A). They predicted steady state radial [Ca 2+ ] gradients extending the full axial distance between T and SR membranes (Fig. 4B,C). Radial dependences of [Ca 2+ ] with distance from the T-SR centre (Fig. 4B) were not altered by the changes in D. Ratios between [Ca 2+ ] at the edge and halfway from the centre, to [Ca 2+ ] at the T-SR centre whether close to T-tubular or SR or between the two membranes, remained at 79% and 14% respectively. Distances for [Ca 2+ ] to fall to 50% between centre and edge values similarly remained unchanged. However, diffusion coefficient value markedly influenced [Ca 2+ ] at the centre of the T-SR junction.
Nevertheless   Table S5). They fall from w = 12 nm to 6.6 nm in hypertonic extracellular solutions 58 . They increase to 20.15 nm and 29.60 nm with fatiguing low-frequency intermittent stimulation 59 and exposure to hypotonic solutions 60 respectively. Computational solutions modelling these variations in w at constant J influx , D and d (Fig. 5) showed early instabilities with time at the greatest w values (Fig. 5A). Nevertheless, all solutions ultimately converged to steady state Ca 2+ microdomains with increased [Ca 2+ ] at the T-SR centre declining with radial distance. However, differences between [Ca 2+ ] values close to the SR, the T-tubular membranes, and within the intervening space occurred at the greater T-SR distances (Fig. 5B). The colourmaps then correspondingly showed marked axial, in addition to radial [Ca 2+ ] non-uniformities (Fig. 5C) with a plume-like tapering. This contrasts with the small axial nonuniformities at w = 12 nm becoming even smaller with its reduction to w = 6.6 nm. This also directly contrasts with the previous near-uniform Ca 2+ microdomain radial profiles through the entire T-SR junctional distance observed in the computations varying J influx and D.
Quantification of these effects as w increased from 6 to 30 nm, showed that close to the T-tubular membrane, [

Ca 2+ microdomains at decreased T-SR diameters.
Thirdly, significant variations in effective areas of membrane appositions occur not only between skeletal muscle T-tubular and SR membrane but also occur in and between other cell types. We quantified these effects successively reducing T-SR junction diameters, d, from the initial 220 nm down to 40 nm, at constant J influx , D and w (Fig. 6, Supplementary Table S6). The computational solutions showed some initial instabilities but ultimately converged even at the smallest T-SR diameters (Fig. 6A). They similarly generated Ca 2+ microdomains in which [Ca 2+ ] declined with radial distance from the T-SR junction centre. These radial gradients accompanied significant axial [Ca 2+ ] variations at the T-tubular and SR membranes and the intervening space (Fig. 6B) resulting in plume-like tapering at the smallest T-SR diameters in the colourmaps (Fig. 6C) (Fig. 7A) indicated the existence of Ca 2+ microdomains with radial concentration profiles (Fig. 7B,C), ratios of [Ca 2+ ] at the edge and halfway between edge and centre, and the centre (79% and 14%), and X' centre , X' 50% , and X' rim values ( www.nature.com/scientificreports/   Fig. 3; Table 1) and resting muscle (Fig. 7) were used to derive the respective proportions of T-tubular membrane facing the T-SR junction model exposed to different microdomain [Ca 2+ ] (Fig. 8). In activated muscle, successively greater proportions of such T-tubular membrane became exposed to successively higher, 0.1 to 10 µM [Ca 2+ ] with increasing depolarization (Fig. 8A)

Discussion
Recent experimental reports implicated hypothetical Ca 2+ microdomains in paradoxical Ca 2+ -mediated effects on skeletal muscle Na v 1.4 activation following pharmacological manoeuvres increasing or decreasing RyR-mediated SR Ca 2+ release 26,27,30 . They went on to suggest that T-SR junctional sites of close membrane proximity, key to excitation contraction coupling, might form diffusion-restricted, ultrastructurally dispersed intracellular subcompartments. Although accounting for only a small proportion of total cytosolic volume, these might sustain regulatory local [Ca 2+ ] heterogeneities in the vicinity of the tubular Na v 1.4 [35][36][37] . These could potentially cause local ~ μM-[Ca 2+ ] alterations previously reported to modify Na + channel function 23 arising from direct Ca 2+ , or indirect, Ca 2+ -calmodulin mediated, binding to Na v 1.4 [18][19][20]23,24 . Our present modelling studies accordingly explored and characterized conditions required for such Ca 2+ microdomain formation within these T-SR junctional structures.
The structural parameters describing sarcomere, surface, T-tubular and SR membrane structure, and distributions, densities and electron microscope ultrastructure of their T-SR junctional regions required for such modelling were available for amphibian skeletal muscle [11][12][13][39][40][41][42][43][44] . [Ca 2+ ] gradients through the resulting formalized geometrical model of a typical T-SR junction in both resting and stimulated muscle fibres were then resolved by finite element method (FEM) solutions of basic Fick diffusion equations. Their boundary conditions first included a RyR-mediated Ca 2+ release into the T-SR space by a uniform Ca 2+ influx density J influx across its SR face. Subsequent Ca 2+ diffusion with a diffusion coefficient established from previous experimental reports from amphibian myoplasm was then mapped through the radially symmetric T-SR junctional space. The second, efflux, boundary condition at the edge of the modelled junction was similarly described by a first order [Ca 2+ ]-dependent process into a well-stirred cytosolic space of infinitely large volume. The latter formalism further matched previous reported eventual first order steady state SERCA-mediated resequestration of the released cytosolic Ca 2+49 .
The boundary conditions simulated conditions both of full and of graded activation by previously reported experimental voltage clamp steps from resting to both 0 mV and varying test potentials and the resulting alterations in cytosolic [Ca 2+ ] 45 . At the influx boundary, Ca 2+ influx densities J influx for each voltage were determined directly from the corresponding reported maximum rates of SR Ca 2+ release, d[Ca 2+ ]/dt, and the model geometrical properties. At the efflux boundary, the resulting [Ca 2+ ] at the edge of the T-SR junction, [Ca 2+ ] edge , was first matched to the corresponding experimental maximum cytosolic concentration [Ca 2+ ] max , by optimising the single free parameter giving exit length λ = 9.2 nm under conditions of full activation by a test step to 0 mV. This assumed the quantities were proportional and close to each other. Both the latter approximations were then further tested in subsequent explorations of varying J influx through different test voltages. In all events, further corrections for any discrepancies arising from a small [Ca 2+ ] edge > [Ca 2+ ] max would tend to enhance rather than reduce the computed [Ca 2+ ] magnitudes.
The resulting T-SR model predicted Ca 2+ microdomains fulfilling the previous suggestions 30 . The computational solutions following step impositions of J influx converged to give steady state T-SR junctional Ca 2+ microdomains. Their heatmap representations demonstrated radial [Ca 2+ ] gradients extending the entire axial T-SR junction distance. The microdomains were quantified radially at the centre, at distances halfway to, and at the edge of the T-SR junction, in axial regions close to the T-tubular and SR membranes and in the intervening T-SR junction space. Furthermore, varying J influx to reflect previously experimentally reported d[Ca 2+ ]/dt, obtained at varying test voltages, at a constant λ value, gave predicted voltage dependences of [Ca 2+ ] edge closely approximating those of the corresponding experimental [Ca 2+ ] max . The accordingly linear, [Ca 2+ ] edge -[Ca 2+ ] max relationship had unity gradient and zero intercept. Furthermore, the [Ca 2+ ] edge values themselves depended linearly on the corresponding Ca 2+ influx Φ influx terms with a gradient fulfilling predictions from the geometrical terms in the equation for the corresponding Ca 2+ efflux.
All these J influx conditions consistently generated T-SR junction Ca 2+ microdomains containing physiologically significant, ~ μM-[Ca 2+ ], heterogeneities graded with imposed depolarization. These could locally elevate [Ca 2+ ] from ~ 0.3-0.4 µM at activation threshold, to 17-20 µM at full activation, relative to the remaining bulk Figure 8. Proportions of T-SR junction T-tubular membrane area exposed to varied tested microdomain [Ca 2+ ]. Results shown for (A) active muscle at different test voltages (Fig. 3) and (B) resting (Fig. 7)  Detailed characteristics of such membrane appositions vary amongst muscle or cell types, and with physiological and physical conditions. Nevertheless, microdomain formation and characteristics were robust through systematic tests at constant J influx that varied Ca 2+ diffusion coefficient, D, T-SR distance, w, and T-SR diameter, d, in turn, holding the remaining variables constant. These tests further made it possible to survey the relative importance of diffusional or geometric properties to microdomain formation and characteristics.
First, alterations in D within the T-SR space could reflect Ca 2+ buffering capacities, κ = (Δ[bound Ca]/Δ[free Ca])), of its contained immobile and mobile buffers, and the diffusion coefficient D mobile of the mobile buffer. Assuming the rapid buffer approximation, the resulting steady state D is related to the free diffusion coefficient D Ca2+ by the expression [65][66][67] : Immobile buffer could reflect fixed Ca 2+ binding sites including negatively charged membrane bilayer phospholipid groups 68 and Ca 2+ -binding domains in Ca 2+ dependent ion channel, cytoskeletal, transport motor and membrane-associated Ca 2+ binding kinase proteins 69 . This would generate local, steady state equilibria between Ca 2+ binding and Ca 2+ diffusion: depleted Ca 2+ -free immobile buffer cannot be replaced by diffusion from remote sites. Immobile buffer would then leave steady-state Ca 2+ microdomains unaffected 67 . Whilst its action could alter the [Ca 2+ ] kinetics, our modelled step impositions of Ca 2+ influxes increased T-SR free [Ca 2+ ] to steady state values over ~ 0.2 ms exponential timecourses. These were 1-2 orders of magnitude shorter than those of experimentally observed Ca 2+ transients 45,70 . In contrast, mobile buffer could influence D to extents dependent upon κ mobile and D mobile . Our computations explored a wide range of conditions extending from limiting maximal D values at infinite Ca 2+ dilution without buffer 55  Secondly, whilst averaging w = 12 nm in resting muscle, axial T-SR distances range from 6 nm with increased extracellular tonicity 58 to 20 nm with fatiguing stimulation 59,60 . Varying w through this range here disrupted microdomain characteristics resulting in heatmaps showing tapering plume-like appearances at the largest T-SR distances. There were increased radial [Ca 2+ ] nonuniformities themselves varying along the axial distance between SR and T-tubular membranes, to extents increasing with increasing w. Nevertheless, [Ca 2+ ] edge varied with w, through an inverse relationship with proportionality constant matching predictions of the T-SR junction model. Thirdly, successive reductions of T-SR junctional diameters, from d = 220 nm to d = 40 nm, similarly disrupted Ca 2+ microdomain heatmaps again giving tapering plume-like forms at the smallest diameters. These were quantified as increased radial non-uniformities and marked axial [Ca 2+ ] differences between regions close to the SR and T-tubular membranes and the intervening space along the axial T-SR distance. Falls in [Ca 2+ ] with radial distance from the T-SR junction centre and [Ca 2+ ] at the SR relative to the T-tubular membranes became less marked with decreasing d. Finally [Ca 2+ ] edge increased with d as expected from the T-SR junction model. Nevertheless, through both these latter modifications in T-SR junction geometry, despite their altered spatial characteristics, the ~ μM-[Ca 2+ ] heterogeneities between their periphery and centre persisted even with more than 100% increases in T-SR distance or 75% reductions in T-SR diameter from control values derived from established morphometric data.
The previous experimental reports had also invoked background, RyR-mediated influxes of Ca 2+ in Ca 2+ microdomain generation in resting in addition to activated muscle 27,30  The present findings taken together could be used to reconstruct the proportions of T-tubular membrane area and therefore of resident Na v 1.4 exposed to successively greater levels of T-SR junction microdomain [Ca 2+ ] in both activated and resting muscle. Successively greater proportions of activated T-tubular membrane became exposed to successively higher, 0.1 to 10 µM [Ca 2+ ] with increasing depolarization. In addition, significant proportions of even resting T-tubular membrane remained exposed to significant, ~ 0.5 µM [Ca 2+ ]. These findings therefore provide a physical basis for the previous suggestions implicating Ca 2+ microdomain formation in observed modifications in Na v 1.4 function 27,30 . Ca 2+ -CaM binding takes place with ~ µM [Ca 2+ ] dissociation constants 72,73 . Feedback µM-[Ca 2+ ] levels arising from SR Ca 2+ release could therefore potentially modify both skeletal 26,27,30 and cardiac muscle 28,29 Na v 1.4 or Na v 1.5 through direct or indirect, Ca 2+ -calmodulin (Ca 2+ -CaM) mediated, actions on their C-terminal domains 16,[18][19][20]23,24 . Such concentrations further match the photoreleaseinduced 1-2 µM cytosolic [Ca 2+ ] elevations previously reported to modify Na v 1.4 function 23 .
In skeletal muscle, elevated T-SR junctional microdomain [Ca 2+ ] could inhibit tubular Na v 1.4 function following normal sustained activity 59,60 and contribute to particular clinical skeletal myopathies 74 . A myotonic hyperexcitability disorder disrupting Ca 2+ -mediated inhibition of Na v 1.4 function has been associated with Na v 1.4 C-terminal EF hand-like domain mutations 75,76 . A condition associated with increased myotube diameters and resting [Ca 2+ ] i , and decreased RyR1-mediated Ca 2+ release reflecting possible abnormalities in triad junction formation and maintenance is associated with another, junctophilin (JP2), mutation 77 . In murine cardiac muscle, Na v 1.5 inhibition followed increased SR Ca 2+ release following pharmacological challenge 28,29 and in RyR2-P2328S variants modelling the pro-arrhythmogenic catecholaminergic polymorphic ventricular tachycardia 78-80 . In these examples, the underlying in vivo source of microdomain Ca 2+ would likely remain the RyR-mediated Ca 2+ release modelled here, as opposed to T-tubular Ca v 1.1 or Ca v 1.2 channel mediated entry of extracellular Ca 2+ . Thus, early Ca 2+ skeletal muscle voltage clamp currents, I Caf (~ 25 µA/cm 2 ) 81 and later cell attached patch clamped cardiomyocyte I CaL (~ 10 pA/pF 82 ; would contribute J influx terms (~ 8.64 × 10 -7 and ~ 6.91 × 10 -8 mol/ m 2 /s assuming similar C T /C S and ξ; respectively) one and two orders of magnitude lower than the corresponding RyR-mediated J influx , at comparable test voltages. The larger skeletal muscle late I Ca (80 µA/cm2, giving J influx ~ 2.76 × 10 -6 mol/m 2 /s) evolves over time courses (hundreds of ms) more prolonged than excitation contraction coupling and shows rapid off kinetics on action potential repolarization 83,84 .
Microdomain µM-[ Ca 2+ ] could also regulate other signalling biomolecules. They are involved in a bell-shaped in vitro open probability relationship for single channel RyR activation and inhibition 14 . Here, cardiac and neuronal, RyR2 and RyR3 are more Ca 2+ -sensitive than RyR1 but all are activated over the ~ 1 µM [Ca 2+ ] predicted in the present analysis 85,86 . However, under their respective in vivo physiological [ATP] and [Mg 2+ ] conditions, cardiac 87 but not skeletal muscle 88 RyR activation involves Ca 2+ -induced Ca 2+ release. Skeletal muscle RyR activation instead involves direct allosteric coupling with either T-tubular Ca v 1.1-L-type Ca 2+ channel voltage sensors 15 or possibly other adjacent SR RyRs themselves coupled to such Ca v 1.1 89,90 . ~ µM Ca 2+ -CaM may also exert other cytosolic effects as on glyceraldehyde 3-phosphate dehydrogenase 91 or itself provide local signaling domains 36 .
Closely apposed membranes potentially mediating localized Ca 2+ signalling involving Ca 2+ -dependent proteins that would similarly permit divergent signalling at different sites also occur in widespread other cell types 9,10 . At smooth muscle SR-plasma membrane appositions 92 , local Ca 2+ could modulate repolarizing Ca 2+ -activated K + channel activity 93 . They also occur in neurons 6 ; cerebellar Purkinje and hippocampal neurones similarly signal using RyR-Ca 2+ release channels [1][2][3][4] . Amongst non-excitable cells, multiple 20-30 nm diameter membrane invaginations in thrombocyte open canalicular systems (OCS) 7,8 form vacuolar structures apposed to membranes of the Ca 2+ -storing deep tubular system (DTS) previously compared with muscle T-SR junctions 94 . These gate inositol trisphosphate receptor rather than RyR mediated Ca 2+ fluxes into the junction upon agonist stimulation. The theoretical analysis here thus combined available experimental anatomical and physiological data and diffusion theory to predict significant [Ca 2+ ] heterogeneities in the skeletal muscle T-SR-junction. Its findings might next prompt investigations of structure and function in these further examples. These might be complemented by direct experimental fluorescent Ca 2+ indicator [Ca 2+ ] TSR measurements [35][36][37] were these to be able to address the small dimensions and dispersed nature of the microcompartments concerned.

Materials and methods
Finite element analysis. Diffusional mechanics and its consequent temporal and spatial [Ca 2+ ] patterns were computationally examined in a model T-SR junction. This involved replicating the geometry within and through which the diffusional processes occurred, pre-processing through meshing and definition of physical conditions including loads, initial and boundary conditions, generation of system solutions and post-processing of the results (Supplementary Fig. S1). The matrix-based programming platform and language MATLAB (version R2020b win64 9.9.0.1467703, version released August 26th, 2020, MathWorks, Cambridge, UK) performed the required data array manipulations and generated all the graphics in Figs. 1, 2, 3, 4, 5, 6, 7, and 8 (https:// www. mathw orks. com/ disco very/ what-is-matlab. html). It was implemented on an IBM compatible computer (CPU: Intel Core; i7-4790 K CPU: 4.40 GHz (4 cores); GPU: ASUS STRIX GeForce GTX970; installed RAM: 16 GB, running Windows (Microsoft, Washington, USA) 10, Home 64-bit version 1909).
The underlying T-SR junction geometry was reconstructed virtually for use in a finite element analysis solving partial differential equations (PDEs) for the resulting diffusional processes with their accompanying boundary conditions (BCs) (See Supplementary File for software archive). The finite element method (FEM) described the complex geometry as a collection of subdomains (elements) by superimposing upon this geometry a mesh of tetrahedral elements joined at their vertices (nodes) and edges. The subdivision accurately represents this complex www.nature.com/scientificreports/ geometry, permits inclusion of dissimilar material properties, and provides a straightforward representation of the total solution whilst capturing local, microdomain, effects.

Boundary conditions and equation solutions. Specified BCs provided values of the field and related
variables, in the present case, the normal derivatives of the field variable in the form of a Neumann-type BC. The FEM equations are formulated such that at the nodal connections, the value of the field variable at any connection is the same for each element connected to the node. Solutions at the edges of each adjacent element are therefore equal, ensuring continuity of field variables between elements, avoiding physically unacceptable gaps or voids in the solution 38 . The original PDE problem is accordingly represented within each element with simpler equations approximating the solution to the original equations. Stationary linear problems whose coefficients are independent of the solution or its gradient yield a linear system of equations. In this case our PDE is time-dependent and hence the system of simpler equations is a set of ordinary differential equations (ODEs) then passed onto MATLAB solvers for numerical integration for solution. The FEM approximates the solution by minimizing the associated error function, automatically finding the linear combination of basis functions closest to the solution u. The FEM could therefore capture both concentration differences local to the T-tubular and SR membranes, and across the entire modelled geometry.
Partial differential equation toolbox™. Meshing and application of the FEM used Partial Differential Equation Toolbox (version 3.5 installed on 8th October 2020 by MathWorks) within MATLAB. This provides functions for solving structural mechanics, heat transfer and general PDEs using the FEM (https:// www. mathw orks. com/ produ cts/ pde. html). PDE Toolbox also provides the ability to automatically mesh the T-SR junction geometry, providing a basis for solving the diffusion PDE, and stores the solution as matrices amenable to various methods of presentation and post-processing of data within MATLAB. The PDE toolbox is designed to solve equations of the form: with a generalised Neumann boundary condition of: where the coefficients m, b, h, f and g can be functions of spatial position, the solution u, or its spatial gradient. In a diffusive system, this generalised problem reduces to the first order equation: with b set to unity, where c represents the diffusion coefficient D, and h represents the boundary flux term for the Neumann condition (compare Eq. 4).
Ethical approval. This entirely theoretical study did not involve animal procedures.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request and are furthermore summarized in the Supplementary file. www.nature.com/scientificreports/