The characteristics of intracellular Ca2+ signals in vivo necessitate a new model for salivary fluid secretion

Salivary fluid secretion involves an intricate choreography to result in the trans-epithelial movement of NaCl and water into the acinus lumen. Current models are based on experimental observations in enzymatically isolated cells where the Ca2+ signal invariably propagates globally and thus appears ideally suited to activate spatially separated Cl and K channels. We monitored Ca2+ signals and salivary secretion in live mice expressing GCamp6F, following stimulation of the nerves innervating the submandibular gland. Consistent with in vitro studies, Ca2+ signals were initiated in the apical endoplasmic reticulum. In marked contrast to in vitro data, highly localized trains of Ca2+ transients that failed to propagate from the apical region were observed. Following stimuli optimum for secretion, large apical-basal gradients were elicited. Given this incompatibility to the previous model, a new mathematical model was constructed to explain how salivary secretion can be efficiently stimulated by apically localized Ca2+ signals.

characteristics of Ca 2+ signals in salivary glands is unknown. Additionally, since only indirect measurements of fluid secretion inferred from volume changes can be made in isolated acinar cell preparations [45], it is unclear how the concentration-dependent characteristics of stimulated Ca 2+ signals in vitro relate to the physiological stimulation of fluid secretion in vivo.
In this study we generated mice expressing a genetically encoded Ca 2+ indicator specifically in exocrine acinar cells and imaged the Ca 2+ signals in submandibular glands in living mice by intravital multiphoton (MP) microscopy. Notably, stimulation of the endogenous neural input to the gland induces Ca 2+ signals with major spatiotemporal properties which are distinct from those observed in isolated acinar cells. Specifically, following moderate stimulation, apically localized Ca 2+ signals are generated which did not propagate globally. Furthermore, the apical signals consist of rapid oscillations. At stimulation intensities optimum for secretion, a significant apical-basal standing Ca 2+ gradient was also established. These striking disparities necessitate a modification of our previous model [34,42] to incorporate the physiological characteristics of Ca 2+ signals in vivo and thus provide a revised framework to better understand salivary fluid secretion.

Results
To generate mice expressing a fluorescent Ca 2+ indicator in exocrine acinar cells, homozygous mice expressing the fast genetically encoded Ca 2+ indicator GCaMP6f [13] floxed by a STOP cassette, were crossed with heterozygous tamoxifen-inducible Mist1 Cre mice [26] . The STOP codon was subsequently excised following tamoxifen gavage (Fig. 1A). The expression of GCaMP6f was investigated in vivo in anaesthetized mice following surgery to expose the submandibular salivary glands (SMG). One SMG was lifted and placed on a small platform with a coverslip placed on top of the gland (Fig 1B). Bipolar stimulation electrodes were inserted to the duct bundle connecting to the imaged SMG tissue, in order to stimulate nerves innervating the SMG. A cover glass was placed on top of the SMG to gently hold the tissue and to maintain a saline solution bathing the gland. MP imaging in serial z sections through the gland revealed expression of GCaMP6f exclusively in acinar cells. GCaMP6f expression was uniformly observed in essentially all acinar cells while GCaMP6f expression was not evident in any other cell type such as blood vessels, nerves, myoepithelial cells or ductal cells. Voids with lack of expression of GCaMP6f were evident in z planes distal from the surface of the gland likely representing the localization of ductal structures (Fig 1C).
We next monitored changes in GCaMP6f fluorescence by intravital microscopy in live animals using MP excitation. Images were acquired at 10 frames per second as described in Methods (each image a 3 frame average). Under basal conditions, prior to stimulation, the cells were quiescent (Fig. 2B). Electrical stimulation at various stimulus strengths delivered to nerves innervating the SMG (1-100 Hz, 12 s initiated at 10s) induced increases in intracellular [Ca 2+ ] in acinar cells ( Fig. 2A,B). Following 1 Hz stimulation, a minority of the cells exhibited small, rapid changes in fluorescence, (n = 4 fields from 4 animals, Fig. 2A,B), while most cells were refractory.
Following higher frequency stimulation (3)(4)(5)(6)(7)(8)(9)(10), the majority of the cells in the field responded and the amplitude of the [Ca 2+ ] signal increased as the stimulation strength was increased ( Fig.   2A,B). During the 12 s of stimulation, overall fluorescence gradually increased and immediately returned to baseline after the termination of the stimulation ( Fig. 2A,B). The maximum [Ca 2+ ] signal evoked during the 12 s stimulation progressively increased and was maintained from 1 to 10 Hz throughout the period of stimulation, but this trend did not hold following higher frequency stimulation (n = 4 fields from 4 animals, Fig. 2C). Following stimulation at 30-100 Hz, the [Ca 2+ ] increase was not maintained for the duration of stimulation. Under these conditions, the [Ca 2+ ] rapidly peaked in the initial 1-3 s then decreased thereafter ( Fig. 2A,B), probably as a result of exceeding the maximal firing rate of the neurons innervating the SMG. Of note, the stimulation likely led to increased neural activity specifically, rather than stimulating acinar cells directly, as a consequence of current leak through the saline depolarizing nerve endings, as the contralateral SMG were unresponsive to even the highest 100 Hz nerve bundle stimulation of the ipsilateral gland (p > 0.1 compared to without stimulation, n=4 fields from 4 animals, Fig. 2C).
To further confirm the physiological consequence of nerve stimulation, we monitored saliva secretion under identical stimulus conditions. There was no significant secretion activated by 1-3 Hz nerve stimulation for 1 min (p > 0.05 compared to the absence of nerve stimulation, n = 4 animals), but saliva secretion progressively increased following 5-10 Hz nerve stimulation (p < 0.05 compared to the absence of stimulation, n = 4 animals) (Fig. 3A). Higher frequency stimulation did not result in augmented secretion, likely reflecting failure to sustain the rise in [Ca 2+ ] ( Fig. 2A). At moderate stimulation intensity (0-10 Hz), there was a linear correlation between the peak [Ca 2+ ] increases by stimulation and saliva secretion (r 2 = 0.995, n = 4 animals, Fig. 3B). These stimulus paradigms (1-10 Hz), therefore likely represent the range of physiological stimulation of both evoked Ca 2+ signals and the resulting saliva secretion and thus characterization of these signals, represent the focus of the remainder of the study.
Following nerve stimulation at low stimulation strengths (1-5 Hz), we observed both the recruitment of responding cells and an increase in the amplitude of the [Ca 2+ ] response (Fig. 4).
To determine more closely how the increased stimulation intensity impacts the characteristics of Ca 2+ responses to result in increased saliva secretion, a software-based work-flow was designed to post-process the image series. Imaged fields were sub-divided into 8 by 8 equal quadrants to yield 64 grid squares (Fig. 4C,D). Each grid was 32 x 32 µm, and typically included a portion of an acinus consisting of 1-5 cells and was used to approximate the spatial and temporal behavior of individual cells or lobules in the field. In order to analyze the activity in each region, the image was then processed to show percent DF/F for each time-lapse frame in the sub-divided regions ( Fig. 4D,E). Using this general scheme, we compared information obtained from 8 x 8 grids, 16 x 16 grids and randomly manually selected regions of interest representing either acinar clusters or single acinar cells. Each scheme produced essentially equivalent data, ultimately validating the use of the 8 x 8 grid (Supplemental Fig. 1). At 1 Hz stimulation, the activity was sporadic, with the responding regions tending to give multiple responses while non-responding regions remained quiescent throughout the stimulus period (Fig. 4F). Overall, 34.0 ± 5.1% (n = 19 fields from 8 animals) of the sub-regions showed significant responses during the stimulation (Fig. 4G).
Moreover, the increase in [Ca 2+ ] in the responding regions exhibited a prolonged latency prior to the initial rise following stimulation (6.72 ± 0.79 s, n = 19, Fig. 4H). Stimulation at higher frequencies progressively recruited more grid regions and shortened the latency to the initial response (Fig. 4G, H). At 10 Hz, virtually all regions responded to the stimulation (98.5 ± 0.6%, n = 19) with minimal latency (1.70 ± 0.27 s, n = 19) (Fig. 4G, H). These data suggest that the increase in overall [Ca 2+ ] and secretion by rising stimulation strength was at least in part due to the recruitment of responding cells and a faster time for cells to respond. Further, the amplitude of the maximal increase in [Ca 2+ ] in the subdivided grids also increased as the strength of the stimulation was augmented (Fig. 4I), indicating that enhancement of individual cellular [Ca 2+ ] responses also contribute to the overall increase of [Ca 2+ ] amplitudes.
We observed that increased stimulus strength resulted in both a higher number of responding regions and a higher peak [Ca 2+ ] amplitude in individual grids. We next examined whether grids which responded at lower stimulus intensity also strongly responded to higher frequency stimulation. The alternative being that the spatial pattern of responding grid regions differed with increasing strengths of stimulus. Grids in an imaging field were ranked based on the peak [Ca 2+ ] evoked by 5 Hz stimulation and then re-sorted based on progressively increasing amplitude ( Fig. 5A-B). We subsequently sorted the data obtained from 1, 3, and 10 Hz stimulation from the identical field in the grid rank order established for 5 Hz stimulation (Fig. 5C). This analysis demonstrated a clear pattern whereby, at each stimulus strength, particular regions of the field of view were most susceptible to stimulation, such that grids giving the largest amplitude response at 1 Hz similarly responded by displaying the largest amplitude signals at 3, 5 and 10 Hz stimulation. (n= 5 fields from 4 animals) (Fig. 5D).
The extent of the increase in [Ca 2+ ] following greater strengths of stimulation was not uniform throughout individual grids. High-responding grid regions showed linearly enhanced increases following stimulation (R 2 = 0.9197, n = 4 fields from 4 animals), while the increase in Ca 2+ signal following increased stimulation in less responsive regions was proportionately smaller While the previous grid-based analysis yielded important information relating to the general properties of the evoked signals within an imaged field, this approach was not suited to provide spatial information on a sub-acinus or lobule scale. To interrogate the sub-cellular spatial characteristics of the Ca 2+ signals evoked by neural stimulation, we generated standard deviation (SD) images from the stimulated portion of the image series. This analysis revealed that at all physiological stimulus strengths the increase in [Ca 2+ ]i was highly heterogeneous, with the largest increases in signal occurring in the apical aspects of acinar clusters (Supp. Fig 2). Figure 7A illustrates a field in which an acinar cluster is highlighted prior to neural stimulation at 5 Hz (red box and image series in Fig 7B). To further analyze the neurally evoked apical Ca 2+ signals, Python software scripts were written in the Jupyter lab notebook environment to automate the designation of apical ROIs from masks generated from average or SD images. These apical ROIs were then subsequently used to quantitate the magnitude and frequency of signals within these regions (see Methods for details). Figure 9A shows an imaging field depicting the stimulated-basal average image for the field generated during the 12s of stimulation at 5 Hz from which 79 apical ROIs were generated from masks shown in Figure 9B/C. Figure 9D (Fig 9F)). In contrast to the robust correlation between saliva secretion and the magnitude of the evoked Ca 2+ signal, the relationship between secretion and the frequency of Ca 2+ oscillation was less strong (R 2 =0.824, Fig 9G).
Although GCaMP6f and fura-2 have similar affinities [7,25] , a possibility exists that the marked differences in the spatiotemporal properties of Ca 2+ signals observed in vivo vs.
documented in vitro are a function of the genetically encoded indicator used in the present studies. We therefore prepared small acinar clumps from Mist1 +/x GCaMP6f +/animals by enzymatic digestion and monitored Ca 2+ signals stimulated by the muscarinic agonist carbachol (CCh) by wide field imaging. At low concentrations of CCh, repetitive Ca 2+ oscillations were evoked ( Fig 10A) with a period of ~0.15 Hz (Fig 10B, n=32 cells, 3 animals). In response to increasing concentrations of CCh, the initial peak tended to increase (Fig 10C, n=32 cells, 3 animals) and the oscillations dampened to generate an elevated plateau (Fig 10A). A statistically significant change in the initial peak height was not evident, possibly reflecting the relatively narrow CCh concentration range tested, or alternatively the relatively high affinity of GCaMP6f.
Oscillations evoked by low concentrations of CCh were initiated in the apical region of the cell and invariably globalized (Fig 10D/E). In total, the spatiotemporal profile of the Ca 2+ responses evoked in acutely isolated GCaMP6f expressing acini were strikingly similar to numerous reports documenting agonist-evoked Ca 2+ responses with conventional chemical Ca 2+ indicators and thus the properties of the genetically encoded indicator per se do not explain the differences observed.
We have previously constructed a series of mathematical models to help understand and explain our experimental results and those from other groups [1,34,42,43,50,51]. Briefly, these multi scale models were based on cyclical global increases in [Ca 2+ ] activating spatially separated ion channels required for sustained secretion in a 3D collection of seven acinar cells forming an acinus. Apical initiation of the Ca 2+ signal results in activation of Clflux through TMEM16a Cl channels and the Ca 2+ signal is subsequently transmitted across the cell cytoplasm by a process of Ca 2+ -induced Ca 2+ release (via RyR or IP3R). This leads to periodic increases of [Ca 2+ ] in the basal aspects of the cell and hence activation of basal Ca 2+ activated K + (KCa) channels. Activation of the KCa channels is critical for maintenance of the electrochemical potential, driving the efflux of Cland thus secretion. The spatiotemporal information from the current in vivo data is obviously not compatible with these models. We have therefore constructed a new model, as shown in

Discussion
An increase in [Ca 2+ ] following neurotransmitter release from parasympathetic neurons is fundamentally important for the stimulation of salivary fluid secretion. It is also recognized that the spatiotemporal properties of intracellular Ca 2+ signals in exocrine cells are pivotal for appropriately activating ion channels necessary for the underlying process [2,20,21]. Despite a wealth of information documenting agonist-stimulated Ca 2+ signals, recorded in dissociated single cells, isolated acini or excised lobules, the properties of physiological Ca 2+ signals in vivo following endogenous neural stimulation have not been reported. We therefore generated transgenic mice that express a genetically encoded Ca 2+ indicator specifically in exocrine acinar cells and used MP microscopy in live mice to document the Ca 2+ signals generated following stimulation of the submandibular nerve which innervates the SMG. By establishing stimulation parameters optimum for fluid secretion, we have inferred the spatiotemporal properties of Ca 2+ signals that ultimately underlie the physiological stimulation of fluid secretion at the lobule level and subsequently at subcellular resolution.
Following threshold stimulation, a small proportion of the imaged field responded after a prolonged latency by eliciting a single, or a few brief transients. These Ca 2+ spikes were initiated in the extreme apical aspects of the cell, consistent with the localization of the majority of IP3R and previous in vitro data [24,34,54]. In stark contrast however, these signals did not globalize and were confined to the apical domain of individual single cells within an acinus. This threshold stimulation resulted in minor measurable fluid secretion. At greater stimulus strengths, where secretion was readily evident, an increasingly larger proportion of the field was recruited with shorter latency to yield multiple Ca 2+ transients with increasing peak magnitude. The increase reflecting an increase in all cells within individual responding acini. Strikingly, these signals were again predominately localized to regions of the cell immediately juxtaposed to the apical PM with limited propagation towards the basal aspects of the cell. Upon optimal stimulation for secretion, essentially the entire field was recruited to produce a sustained increase with periodic fluctuations superimposed on the plateau after a minimal latent period. Notably, while the affinity of the indicator may have masked a greater magnitude of apical Ca 2+ rise, or more prominent fluctuations in the signal at these higher stimulus strengths, a prominent apical to basal gradient was readily evident. Further increases in the intensity of the stimulation resulted in a rapid peak followed by a gradual waning of the signal and consistent with the absence of a sustained Ca 2+ signal, the rate of fluid secretion was diminished. This reduced secretion likely results from a "tetanic" stimulus where either neurons fail to fire during their refractory period, neurotransmitter release is exhausted or acinar cell Ca 2+ signaling is perturbed, during this likely non-physiological stimulation. Over the range of physiological stimulus strengths, the extent of fluid secretion was strongly correlated with the magnitude of the peak Ca 2+ rise, rather than necessarily the frequency of Ca 2+ transients. This conclusion is consistent with our previous modeling and experimental work that has indicated that the integrated increase in [Ca 2+ ] during stimulation is the primary driving factor for efficient fluid secretion [42].
While local, non-propagating Ca 2+ signals are reported in isolated pancreatic acinar cells, notably in whole-cell patch clamped cells where Ca 2+ buffering is imposed [46,47], the spatial characteristics of the Ca 2+ signal observed in vivo in the present study represents a major Conceivably, alterations in the polarization of the ER, such that IP3R or SERCA are not maintained in the extreme apical ER [23,24] , or a change in the distribution of mitochondria [11] to disrupt ER-mitochondrial junctions and effective Ca 2+ sequestration from release sites could contribute to the altered Ca 2+ signals observed in vitro.
A further important consideration when interpreting the characteristics of these Ca 2+ signals is the obvious difference in delivery of the secretagogue in these experiments. In contrast to superfusion resulting in a defined "square" pulse of stimulating agent, as commonly used in in vitro studies, where in the acinar cells likely experience a complex, stimulus intensity, and timedependent gradient of neurotransmitter. In our studies, the concentration of ACh released from neurons would be predicted to rapidly increase and then decay as it is hydrolyzed by acetylcholinesterases. As the frequency of stimulation is increased, the residual concentration of ACh present in the vicinity of acinar cells would be predicted to increase with a "saw-tooth" profile as the amount of ACh release exceeds the rate of degradation following more frequent, repetitive cycles of neurotransmitter exocytosis. In this scenario, when a threshold concentration of neurotransmitter is reached that generates sufficient IP3, an increase [Ca 2+ ]i, will be triggered.
The time to reach this threshold would define the latent period for that particular cell and the relative sensitivity to neurotransmitter will predicate whether an individual cell is responsive to a given stimulation. At optimum stimulation frequencies the residual ACh concentration would further increase and remain above the threshold for the majority of acinar cells to respond. This would facilitate greater Ca 2+ release manifested as an increased peak and Ca 2+ oscillations by the entire field. A summary of our observations and these ideas is presented in Supp Fig. 3. Notably, neural stimulation does not simply "pace" cells, as the frequency of oscillations does not correlate with the stimulus frequency. It is therefore likely, that the oscillations result from the inherent properties of IP3R and SERCA pumping activity but where the kinetics are influenced by the dynamic changes in neurotransmitter concentrations. This may contribute to the differences in oscillation frequency observed when comparing in vivo and in vitro data.
We observed considerable spatial heterogeneity in the responsiveness and peak amplitude of cells across the field of view. Specifically, there were sections of the field which exhibit the greatest sensitivity and largest responses at all stimulus intensities. Conceptually, this may reflect acini with expression of a higher density of cell surface receptors. Alternatively, highly responsive regions may represent cells physically closer to synapsing neurons and thus experiencing higher concentrations of ACh. The innervation of SMG lobules visualized in a reporter mouse expressing fluorescent protein in parasympathetic neurons has been reported to be fairly homogeneous [41]. Therefore, a further possibility is that spatial homogeneity arises in part because increasing numbers of axons are recruited to fire as stimulus intensity increases thus directly innervating a larger proportion of cells.
Clearly our previous models describing how Ca 2+ signals promote salivary secretion are not consistent with the new in vivo data presented here. Fundamentally, if there is not an obligatory global intracellular Ca 2+ wave, there is no means of activating the KCa channels in the basal PM that maintain salivary secretion. We have, however previously shown that the KCa channels, KCa3.1 ("IK") and KCa1.1 ("BK") together with Na/K ATPases are present in the apical membrane of salivary acinar cells [1,2]. At the time these results were puzzling, as there was no obvious reason why the apical membrane should express either of these transport proteins.
However, in combination with the new data presented here, the physiological importance of the apical KCa and Na/K ATPases is revealed. If [Ca 2+ ] increases only in the apical region of the cell there must be KCa channels in the apical region, in order to maintain the electrochemical driving force for Clsecretion. On the other hand, since the primary saliva has a low [K + ], this K + must be removed from the lumen. This is the proposed function of the Na/K ATPase in the apical membrane.
Taken together, the present studies describe the spatial temporal characteristics of physiological Ca 2+ signals driving salivary fluid secretion. Further, our findings highlight that caution should be exerted in extrapolating conclusions from ex vivo studies to physiological Ca 2+ signals and function in vivo. Finally, it is envisioned that the present studies will provide a framework for investigating if Ca 2+ signaling is disrupted in disease states such as Sjögren's syndrome or radiation-induced xerostomia which result in pronounced hyposecretion of saliva [9,45].

Methods and Materials
Tamoxifen, corn oil, ketamine and Xylazine were purchased from Sigma Chemical, St Louis, Mo.
Physostigmine was purchased from Tocris Minneapolis, Mn.

Generation of mice expressing GCaMP6f in exocrine acinar cells.
Adult female GCaMP6f flox mice (Jackson Laboratory 029626) were crossed with male heterogenous Mist1 CreERT2 (Jackson Laboratory 029228, kindly donated by Dr. Catherine Ovitt, University of Rochester) to generate Mist1 +/x GCaMP6f +/transgenic mice. Tamoxifen (0.25 mg/g body weight) dissolved in corn oil was administered to mice at least 6 weeks old of either gender, by oral gavage once a day for three consecutive days. Imaging was carried out 1-4 weeks after the last tamoxifen administration. All animal procedures were approved by University Committee on Animal Resources (UCAR-2001-214E).

In vivo MP imaging
The animal was anesthetized with Ketamine (80 mg/kg body weight, i.p.) and Xylazine (10 mg/kg body weight, i.p.). The animal was restrained with the ventral side up, and the animal's body temperature was maintained at 37ºC by a heat pad during surgery and imaging. A submandibular gland was exposed by making an incision in the skin on the ventral side of the neck. Connective tissue around the gland was teased away using forceps to allow the gland to be raised away from the body but remaining attached by a duct/nerve/blood vessels bundle. A pair of tungsten wires (WPI, Inc) was inserted to the bundle as stimulation electrodes. The lifted gland was placed on a 10 x 13 mm custom build small holder situated directly above the neck so that the position of the gland was not influenced by movement as a result of breathing. The gland was held by a cover glass to flatten the surface and to keep tissue moist with Hanks salt solution between the holder and the cover glass ( Figure 1B). A 25x water immersion lens (Olympus XLPlan N 1.05 W MP) equipped with an objective heater (OKOLab COL2532) kept at 37ºC was utilized for MP imaging.
An upright 2-photon microscope system (FVMPE-RS, Olympus) with an excitation laser (InSight X3, Spectra-Physics) tuned at 950 nm and emission collected at 495-540 nm was utilized. Images were captured using Olympus FV31s-SW software at 30 fps framerate then averaged every 3 frames to achieve 10 fps image collection rate. Imaging depths between 10 and 55 µm from the surface of the gland were routinely utilized. PMT settings were fixed at 600V, 1x gain, and 3% black level, with excitation laser power adjusted per animal and according to imaging depth.
Stimulation was generated by a stimulus isolator (Iso-flex, A.M.P.I.) set 5 mA, 200 µs, at the indicated frequency with train frequency and duration controlled by a train generator (DG2A, Warner Instruments). Images were captured typically for 30 s continuously. At least 1 min interval was given before imaging in the same field at a different stimulation strengths. In selected experiments, physostigmine (0.1 mg/kg body weight, i.p.) was administered to the animal on the microscope stage, at least 15 min prior to the imaging of post-physostigmine Ca 2+ signals.

Saliva secretion measurement
An anesthetized animal with stimulation electrodes attached to the duct bundle were gently held with ventral side up. A piece of filter paper, roughly 2 x 10 mm, were weighed before placing it in the animal's mouth. A stimulation of 0-100 Hz was given to a gland for 1 min, and the paper was immediately removed and weighed. The difference of the weight of the paper represents the amount of saliva secreted out in the animal's mouth. A new paper was placed in the mouth for each stimulation train.

Acinar cell isolation and imaging
SMG acinar cells were enzymatically isolated from 2-4 month old, Mist1 +/x GCaMP6f +/mice of both sexes. To isolate acinar cells, glands were extracted, connective tissue was removed, and glands were minced. Cells were placed in oxygenated dissociation media at 37˚C for ~30 minutes with shaking. Dissociation media consisted of Hank's Balanced Salt Solution containing CaCl2 and MgCl2 (HBSS), bovine serum albumin (0.5%), and Collagenase Type II (0.2 mg/mL, Worthington).
Cells were washed twice in HBSS with 0.5% BSA and resuspended in a HBSS solution containing 0.5% BSA and 0.02% trypsin inhibitor. Cells were then resuspended in imaging buffer (10 mM HEPES, 1.26 mM Ca 2+ , 137 mM NaCl, 4.7 mM KCl, 5.5 mM glucose, 1 mM Na2HPO4, 0.56 mM MgCl2, at pH 7.4) and seeded onto a coverslip to allow attachment of cells. Cells were then perfused with imaging buffer and stimulated with agonist. Ca 2+ imaging was performed using an inverted epifluorescence Nikon microscope with a 40 X oil immersion objective (NA=1.3). Cells were excited at 488 nm from a monochromator, and emission was monitored at 530 nm. Images were captured every 500 ms with an exposure of 10 ms and 4 × 4 binning using a digital camera driven by TILL Photonics software. Image acquisition was performed using TILLvisION software and data was exported to Microsoft excel.

Batch Image analysis
Amplitude, latency, and percent of responding areas in fluorescence fields were analyzed using ImageJ software (NIH). XY drift during the imaging, when present, was corrected by applying the Descriptor-based series registration plugin (Stephan Preibisch). The imaging field was divided into 8 by 8 grids to yield 64 regions of interest (ROIs) of dimension 32x32 µm square. The average intensity of each ROI grid was generated for each frame. The average of first 100 frames prior to stimulation served as baseline fluorescence (F0), and %DF/F0 ((F-F0)/ F0 x 100) was calculated using the Image calculator function, so that the converted 32-bit image series represents [Ca 2+ ] changes over time expressed as %DF/F0 in 8x8 alley of grids in the XY dimension (Fig. 5A). The standard deviation of the baseline 100 frames from grids in each image series provided an estimate of the noise level. A change greater than 4 times the standard deviation from the F0 value was considered as a stimulated Ca 2+ signal and the ROI therefore designating a responding grid. To calculate latency prior to a response the time at the first incident Ca 2+ signal meeting the criteria as a responding grid in each ROI after the initiation of the stimulation were recorded.
Non-responding grids were excluded from this analysis. Statistical analyses were performed with paired t test, one-way ANOVA, and linear regression using Prism (GraphPad).

Sub cellular image analysis
A computer based, automated region-of-interest (ROI) detection method that specifically targets cellular regions with pronounced changes was developed since no off-the-shelf tool met our particular requirements to efficiently process the data. The new software toolset consists of an assembly of Python (https://www.python.org) scripts deployed in a collection of Jupyter Lab (https://jupyter.org) "notebooks". The notebooks were designed to automate repetitive analysis steps and give flexibility to interactively explore and analyze the image data sets and makes use of image processing algorithms in the scikit-image (https://scikit-image.org) package. ROI were identified by generating a difference image whereby the average image prior to stimulation was subtracted from the average image during the period of stimulation on a pixel-by-pixel basis ( Figure 10A). An initial mask was created by simple pixel intensity thresholding ( Figure 10B). Subsequent filtering by binary dilation and binary erosion removes undesired small regions and smooths the remaining regions ( Figure 10C). The notebook then generates plots of pixel intensity over time for each ROI generated and utilizes the same mask for images sets generated from the same field at differing frequencies of stimulation ( Figure 10D). A further script block was generated to objectively automate the identification of peaks generated by stimulation. For peak processing several signal processing operations from the scipy (https://www.scipy.org) package were employed. The script employs a sequence of signal resampling, low-pass zero phase filtering, high-pass zero phase filtering and a generic peak detector followed by a mapping back into the original response data. Detected peaks for each stimulation frequency were plotted in the notebook as black dots overlaid on region summary plots as shown in Figure 10D.

Model details
The model has two interconnected modules; a Ca 2+ oscillation module and a fluid secretion module. The two modules are coupled via changes in cell volume, and by the activation of KCa and ClCa by Ca 2+ . The Ca 2+ oscillation module is a reaction-diffusion equation and is solved in a three-dimensional domain reconstructed from experimental structural data [51]. A finite element method was used. All other ions are assumed to be homogeneously distributed in the cell, and are thus described by a system of ordinary differential equations. The individual models for the various transporters, channels and exchangers are all given in [51]. The model equations were solved in a single, uncoupled cell (cell 4 in the notation of [51]). Similar results were obtained from all the other cells). The volume of the cell is taken into account, not by a full remeshing at each time step, but simply by a volume scaling factor in the reaction diffusion equation.
The Ca 2+ oscillation module is based on a closed-cell Class I model in which oscillations arise via sequential activation and inactivation of the IP3R. [Ca 2+ ] buffering is assumed to be fast and linear, and the IP3R model is taken from [43]. Ca 2+ obeys a reaction-diffusion equation in the cell interior, but the effective diffusion coefficient of Ca 2+ is small (due to buffering) which allows for localised increases in [Ca 2+ ] which are not propagated throughout the cell. Because the IP3R are situated no more than 50 nm from the ClCa (which are on the apical membrane), release of Ca 2+ through the IP3R is modelled as a boundary term, which avoids the necessity for a highresolution finite element mesh in the apical region, which would greatly slow the calculations.
The diffusion coefficient of IP3 is two orders of magnitude greater than that of Ca 2+ , and thus IP3 is effectively spatially homogeneous throughout the cell.
In the new model, all parameters and equations remain unchanged from those in the old model [51] with the following exceptions.
1. The maximum Ca 2+ flux through RyR (VRyR) was set to 0. This is the simplest way of removing RyR from the model.

2.
A K channel, with the same conductance as the basal KCa, was introduced into the apical membrane, with the equations for cytosolic and lumenal [K + ] altered to compensate. The apical/basal K + current ratio was determined by the apical/basal surface area ratio.
3. Na/KATPases were introduced into the apical membrane, using the same model and parameters as in (ref). The only change is that the total Na/KATPase activity was assumed to be unchanged, with 70% occuring in the basal membrane, and 30% in the apical membrane. The equations for cytosolic and lumenal [K + ] and [Na + ] were altered to compensate. 4. VPLC , which controls the level of agonist stimulation, was set to 0.01 µM/s. As before [34, 50,51], PLC activity was assumed to be dependent on [Ca 2+ ] (although this has little effect on the oscillations), and nonzero only at the basal membrane.           Fig. 3A), which showed a linear regression with R 2 = 0.824 (black line). ** P<0.01 ANOVA with Tukey test.   11. A proposed updated model for salivary secretion. The mathematical model is based on the fluxes and processes shown here. In contrast to previous models, the apical region now contains all the machinery needed for saliva secretion, included KCa channels and Na/K ATPases.
Ca 2+ is released predominately in the apical region, from IP3R that are situated in close proximity to TMEM16a, and there is no requirement for a propagated wave of increased [Ca 2+ ] across the cell. responses in grey are much lower than the apical responses in red.